Suppose that you have a system of PDEs to solve. At least for simplicity, let's assume it's time independent, quasi-linear (linear in its derivatives) solved on a rectangular grid in (x,y) space, and with boundary conditions specified all around. My question is more general, but let's start here.
There might be two dependent variables, $u(x,y)$ and $v(x,y)$. The general equation might of the form:
$$ a(x,y) Y_{xx} + b(x,y) Y_{yy} + cY_{xy} + d(x,y) Y_x + e(x,y) Y_y = f(x, y, Y) $$
where all the functions $a$ through $e$ are 2x2 matrices, $f$ is a 2x1 matrix, and $Y$ is
$$ Y(x,y) = \begin{pmatrix} u(x,y) \\ v(x,y) \end{pmatrix} $$
Suppose that you would like to compute a finite-difference numerical approximation. Assume that the grid-points are evenly spaced in x and y. You might discretize $x$ into $N$ points and $y$ into $M$ points. And then you would construct a solution as a 2(nm) column vector.
$$ X = \begin{pmatrix} u_{11} \\ \vdots \\ u_{n1} \\ \vdots \\ u_{nm} \\ v_{11} \\ \vdots \\ v_{nm} \end{pmatrix} $$
You would then solve some sort of matrix equation $$ (A + B + C + D + E)X = b, $$
where the 2nm x 2nm matrices $A$ through $E$ are finite difference matrices for the corresponding differential operators. The combined matrix in the brackets will have some sort of block banded structure. Their edge values may be complicated due to boundary conditions.
My question is simple: is there an easy way (particularly in Matlab) to generate the matrices $A$ through $E$? Or at the very least, a simple guide to creating all the necessary matrices (except perhaps the modifications you need for boundary conditions, which might be a manual input).
The problem is that I can do this by hand, but it's (i) a lot of algebra; and (ii) subject to a lot of typos and errors when you implement.
For instance, here is a code that gives you all the necessary vectors for a 1D finite difference 'matrix' for different orders of derivatives. Is there an analogy for systems of PDEs on the plane?
