Matlab help related with discretization of second order elliptic partial differential equation

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. • Is this a homework or other academically assigned problem? Commented Apr 1, 2014 at 22:36 1 Answer You will have to augment whatever matrix you get as a result of applying a finite difference stencil discretizing the PDE with some equations representing the boundary conditions. The reason for this can be seen from looking at the one dimensional case. Say you had a just the ODE u''(x) = 0 on the domain x \in (0,1). Say you want to discretize the domain into 5 intervals. 0 --- 0.2 --- 0.4 --- 0.6 --- 0.8 --- 1.0 Above I made a schematic diagram of the discretized domain. At the internal points you can get a finite difference approximation of the derivative as \frac{u_{n-1}-2 u_n + u_{n+1}}{\Delta x} where n is the index of the point the stencil is centered at. So you can get equations for the x=0.2 point where say n=1 and the x= 0.4 point where say n=2 up to the x = 0.8 point where n=4. But for n=0 and n=5 it won't work because the stencil would want to look u_{-1} or u_{6} which are outside our domain. So here is where the boundary conditions come in. But for the internal points you can get 4 equations. In matrix form they would be$$\left( \begin{array}{cccccc} 1 & -2 & 1 & 0 & 0 & 0 \\ 0 & 1 & -2 & 1 & 0 & 0 \\ 0 & 0 & 1 & -2 & 1 & 0 \\ 0 & 0 & 0 & 1 & -2 & 1 \\ \end{array} \right).\left( \begin{array}{c} u_0 \\ u_1 \\ u_2 \\ u_3 \\ u_4 \\ u_5 \\ \end{array} \right)=\left( \begin{array}{c} 0 \\ 0 \\ 0 \\ 0 \\ \end{array} \right)\Delta x$$Adding the boundary conditions would make this a well posed problem. For example if u(0) = 1 and u(1)=5$$\left( \begin{array}{cccccc} 1 & 0 & 0 & 0 & 0 & 0 \\ 1 & -2 & 1 & 0 & 0 & 0 \\ 0 & 1 & -2 & 1 & 0 & 0 \\ 0 & 0 & 1 & -2 & 1 & 0 \\ 0 & 0 & 0 & 1 & -2 & 1 \\ 0 & 0 & 0 & 0 & 0 & 1 \\ \end{array} \right).\left( \begin{array}{c} u_0 \\ u_1 \\ u_2 \\ u_3 \\ u_4 \\ u_5 \\ \end{array} \right)=\left( \begin{array}{c} 1 \\ 0 \\ 0 \\ 0 \\ 0 \\ 5 \\ \end{array} \right)$\$

It gets a little more complicated in two dimensions because you have to index the solutions values with a double index, but still list them in a one dimensional vector. Anyway if you posted an attempt you've made we could help you with it, but it's probaly unreasonable to ask people to write the matlab code for you from scratch.