\(\nabla ^{2}u=f\) in 2D is deﬁned as \(\frac {\partial ^{2}u}{\partial x^{2}}+\frac {\partial ^{2}u}{\partial y^{2}}=f\) and the Laplacian operator using second order standard diﬀerences results in \(\left ( \frac {u_{i-1,j}+u_{i+1,j}+u_{i,j-1}+u_{i,j+1}-4u_{i,j}}{h^{2}}\right ) =f_{i,j}\) where \(h\) is the grid size. The above is solved for \(x\) in the form \(Ax=f\) by generating the \(A\) matrix and taking into consideration the boundary conditions. The follows show how to generate the sparse representation of \(A\). Assuming the number of unknowns \(n=3\) in one direction, there are 9 unknowns to solve for and the \(A\) matrix will be \(9\) by \(9\).