I want to solve, in MatLab, a linear system (corresponding to a PDE system of two equations written in finite difference scheme). The action of the system matrix (corresponding to one of the diffusive terms of the PDE system) reads, symbolically (`u`

is one of the unknown fields, `n`

is the time step, `j`

is the grid point):

and fully:

The above matrix has to be intended as `A`

, where A*U^n+1 = B is the system. `U`

contains the 'u' and the 'v' (second unknown field of the PDE system) alternatively: U = [u_1,v_1,u_2,v_2,...,u_J,v_J]. So far I have been filling this matrix using `spdiags`

and `diag`

in the following expensive way:

```
E=zeros(2*J,1);
E(1:2:2*J) = 1;
E(2:2:2*J) = 0;
Dvec=zeros(2*J,1);
for i=3:2:2*J-3
Dvec(i)=D_11((i+1)/2);
end
for i=4:2:2*J-2
Dvec(i)=D_21(i/2);
end
A = diag(Dvec)*spdiags([-E,-E,2*E,2*E,-E,-E],[-3,-2,-1,0,1,2],2*J,2*J)/(dx^2);`
```

and for the solution

```
[L,U]=lu(A);
y = L\B;
U(:) =U\y;
```

where `B`

is the right hand side vector.

This is obviously unreasonably expensive because it needs to build a JxJ matrix, do a JxJ matrix multiplication, etc.

Then comes my question: is there a way to solve the system without passing MatLab a matrix, e.g., by passing the vector `Dvec`

or alternatively directly `D_11`

and `D_22`

? This would spare me a lot of memory and processing time!