I am reading this paper. In Example 2 from this paper, linear system of equation $Ax = b$ is given, where coefficient matrix $A$ has been generated by he five-point discretization of the following second order elliptic partial differential equation:
$-\frac{\partial}{\partial x}(a\frac{\partial u}{\partial x})-\frac{\partial}{\partial y}(b\frac{\partial u}{\partial x}) + \frac{\partial}{\partial x}(cu)+ \frac{\partial}{\partial y}(du) +fu = 0 $$~~~~ (1)$
with $a(x, y)>0$, $b(x, y)>0$, $c(x, y)$, $d(x, y)$ and $f(x, y)$ defined on a unit square region $\omega = (0, 1)\times (0,1)$, and Dirichlet boundary condition $u(x, y) = 0$ on boundary of $\omega$.
Now for the equation $(1)$ with $a(x, y) = b(x, y) = 1$, $c(x, y) = \cos(x/6)$, $d(x, y) = \sin(y/6)$, $f(x, y) = 1$ with three uniform meshes of $hx = hy = 1/11, 1/21, 1/31$ we get three matrices of orders $100\times 100$, $400\times 400$, $900\times 900$, where $hx$ and $hy$ refer to the mesh sizes in the $x$ direction and $y$-direction, respectively.
I need help to generated matlab code for this problem. I have no idea how to apply matlab to generate these matrices. I would be very much grateful for any kind of help.
Thank you very much for your time.