Efficient implementation of D x / D.T @ v

In several solvers, sparse matrix is used to implement the finite difference operator as
```
len_y = len(self.y)
data = np.array([np.ones(len_y), -np.ones(len_y)])
diags = np.array([0, 1])
D = spdiags(data, diags, len_y-1, len_y)
```
I believe it could be implemented with a forward operator `np.diff` and an adjoint `-np.diff(x, append=0, preprend=0)`
A quick profiling give me x2 - x3 speedup, but most importantly it generalizes better if one wants to do 2D.

No file or test is named in the issue. Locate the sparse finite-difference construction across the solvers, then compare the proposed forward and adjoint operators for numerical equivalence and performance. Done means the relevant solvers use the equivalent efficient implementation without changing their results.

Written by the indexing model from the issue text.