For $\Omega = (0,1)^2 \subseteq \mathbb R^2$ and $f \in C(\Omega)$ consider the boundary value problem:
$- \Delta u(x,y) = f(x,y)~ \forall (x,y) \in \Omega,~ u(x,y)=0 ~ \forall (x,y) \in \partial \Omega$
For discretization let be: $n \in \mathbb N, ~ h= \frac{1}{n+1}, ~ x_i = i \cdot h, ~ y_j = j \cdot h$ for $i,j= 0, \ldots, n+1$.
The task is to find the linear equation system for the discretized solution: $u_{i,j} = u(x_i,y_j), ~ i,j=1, \ldots, n$.
Therefore I have to deduce an approximation of $\Delta u(x_i,y_j)$, dependent on $u_{i,j},u_{i+1,j},u_{i-1,j},u_{i,j-1},u_{i,j+1}$, from taylor expansions of $u_{i-1,j},u_{i+1,j},u_{i,j-1},u_{i,j+1}$ around the point $(x_i, y_j)$ and then use the result to deduce a linear equation system for $(u_{1,1}, \ldots, u_{1,n}, u_{2,1}, \ldots, u_{2,n}, \ldots, \ldots, u_{n,1}, \ldots u_{n,n})$ where the right side is given by discretization $f_{i,j} = f(x_i, y_j)$.
Unfortunately I don't even know how to start here and would be glad for any help.
Thanks in advance!
You're on the right track. \begin{align} \Delta u(x,y) = \partial_{xx} u(x,y) + \partial_{yy} u(x,y) \end{align} Now the most common way is to apply finite differences to the second derivatives, i.e. \begin{align} \partial_{xx} u(x_i,y_j) \approx \frac{1}{h^2}[u(x_{i-1},y_j)-2u(x_i,y_j)+u(x_{i+1},y_j)] \\ \partial_{yy} u(x_i,y_j) \approx \frac{1}{h^2}[u(x_{i},y_{j-1})-2u(x_i,y_j)+u(x_{i},y_{j+1})] \end{align} And your Laplacian becomes \begin{align} \Delta u(x_i,y_j)=\frac{1}{h^2} [u(x_{i-1},y_j)+u(x_{i+1},y_j)-4u(x_i,y_j) +u(x_{i},y_{j-1})+u(x_{i},y_{j+1})] \end{align} Or by writing $u(x_i,y_j)$ as $u_{i,j}$ \begin{align} \Delta u_{i,j}=\frac{1}{h^2} [u_{i-1,j}+u_{i+1,j}-4u_{ij} +u_{i,j-1}+u_{i,j+1}] \end{align} Now do this for all $i,j$ and set up a system of linear equations $Au=f$ where $u$ contains all $(n+1)^2$ entries and $f$ the right handside. Do you know, what $A$ looks like?