How to solve parabolic equation via implicit Euler in 2 dimensions? I have the following parabolic equation: 
$$
\frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2} + \frac{\partial ^2u}{\partial y^2}
$$
over domain  $(x,y)\in [0,10] \times [0,10]$ where $\Delta x = \Delta y = 1 \times 10^{-2}$ and $\Delta t = 2.4 \times 10^{-3}$. Initial and boundary conditions are:
$$
\left\{
    \begin{array}{ll}
    u(x,y,0) = 10\\
    u(0,y,t) = 4\\
    u(10,y,t) = 4\\
    u(x,0,t) = 9\\
    u(x,10,t) = 9\\
    \end{array}
\right.
$$
I tried to rewrite the equation as follows:
$$
\frac{u_{i,j}^{k+1}-u_{i,j}^{k}}{\Delta t} = \frac{u_{i+1,j}^{k+1}-2u_{i,j}^{k+1}+u_{i-1,j}^{k}}{\Delta x^2} +  \frac{u_{i,j+1}^{k+1}-2u_{i,j}^{k+1}+u_{i,j-1}^{k}}{\Delta y^2}
$$
I initially wanted to write the equation as a tridiagonal matrix:
$$
-\frac{1}{\Delta x^2}(u_{i+1,j}^{k+1}+u_{i-1,j}^{k+1}) - \frac{1}{\Delta y^2} (u_{i,j+1}^{k+1}+u_{i,j-1}^{k+1})+(\frac{1}{\Delta t}+\frac{2}{\Delta x^2}+2\frac{1}{\Delta y^2})u_{i,j}^{k+1}= \frac{u_{i,j}^{k}}{\Delta t}.
$$
However, I didn't found a way to do so. I am not sure how to solve this, but is it possible that I have to solve this by iteration like:
$$
u_{i,j}^{k+1}= \Big(\frac{u_{i,j}^{k}}{\Delta t}+\frac{1}{\Delta x^2}(u_{i+1,j}^{k+1}+u_{i-1,j}^{k+1}) + \frac{1}{\Delta y^2} (u_{i,j+1}^{k+1}+u_{i,j-1}^{k+1})\Big)(\frac{1}{\Delta t}+\frac{2}{\Delta x^2}+2\frac{1}{\Delta y^2})^{-1}.
$$
 A: As user_of_math suggested, you should be consistent with your $k$ index.Let $k$ denote the index for time, $i$ for $x$-direction in space and $j$ for $y$-direction in space. Then you use the wrong discretization. It should always be $k+1$ as the time index except for the time derivative since your scheme should be implicit.
\begin{align}
\frac{u^{k+1}_{i,j}-u^{k}_{i,j}}{\Delta t}
=
\frac{u^{k+1}_{i+1,j}-2u^{k+1}_{i,j}+u^{k+1}_{i-1,j}}{\Delta x^2}
+
\frac{u^{k+1}_{i,j+1}-2u^{k+1}_{i,j}+u^{k+1}_{i,j-1}}{\Delta y^2}
\end{align}
or when you want to solve for the next time step
\begin{align}
u^{k+1}_{i,j}
=
u^{k}_{i,j}+\Delta t\Bigl(\frac{u^{k+1}_{i+1,j}-2u^{k+1}_{i,j}+u^{k}_{i-1,j}}{\Delta x^2}
+
\frac{u^{k+1}_{i,j+1}-2u^{k+1}_{i,j}+u^{k+1}_{i,j-1}}{\Delta x^2}\Bigr)
\end{align}
and since $\Delta x =\Delta y$
\begin{align}
&u^{k+1}_{i,j}
=
u^{k}_{i,j}+\frac{\Delta t}{\Delta x^2}\Bigl({u^{k+1}_{i+1,j}+u^{k+1}_{i-1,j}}
-4u^{k+1}_{i,j}
+
{u^{k+1}_{i,j+1}+u^{k+1}_{i,j-1}}\Bigr)
\\
\Leftrightarrow&
u^{k+1}_{i,j}-\frac{\Delta t}{\Delta x^2}\Bigl({u^{k+1}_{i+1,j}+u^{k+1}_{i-1,j}}
-4u^{k+1}_{i,j}
+
{u^{k+1}_{i,j+1}+u^{k+1}_{i,j-1}}\Bigr)
=u^k_{i,j}
\end{align}
You get this for each pair $(i,j)$ giving you $N_x \times N_y$ equations. Assume $N=N_x=N_y$. Now you need to solve a linear system in each time step. You have to come up with some ordering for your matrix and right hand side. The standard ordering would give you a matrix $A \in \mathbb{R}^{N^2\times N^2}$ which looks like this:
\begin{align}
A = 
\begin{pmatrix} 
\tilde A &-\tilde I &0 & ...&...&0 \\
-\tilde I & \tilde A&-\tilde I &0 & .. & 0 \\
0 &-\tilde I & \tilde A &-\tilde I & 0 &..\\
 &&& \ddots \\
0 &...&...&0&-\tilde I & \tilde A
\end{pmatrix}
\end{align}
Where $\tilde I$ is the $n\times n$ identiy matrix, scaled by $c:=\Delta t / \Delta x^2$ and $\tilde A$ is an $n\times n$ matrix.
\begin{align}
\tilde A = 
\begin{pmatrix} 
1+c&-c &0 & ...&...&0 \\
-c  & 1+c&-c  &0 & .. & 0 \\
0 &-c & 1+c &-c  & 0 &..\\
 &&& \ddots \\
0 &...&...&0&-c  & 1+c
\end{pmatrix}
\end{align}
Now you solve 
\begin{align}
Au^{k+1}=u^k
\end{align}
during each time step.
Edit: Of course you have to take care of the boundary values during the computation. I did not consider this.
