I have the following equation:
$$ \nabla^2u = f $$
over $\Omega: [0,10] \times [0,10]$ where boundary conditions:
$$ \left\{ \begin{array}{ll} \frac{\partial u (0,y)}{\partial x} = 0 \\ \frac{\partial u (x,10)}{\partial y} = 1 \\ u(10,y) = 5 \\ u(x,0) = 10 \\ f = 0 \end{array} \right. $$
The stepsize is chosen $\Delta x=\Delta y$. I am not sure whether this equation was solvable by way of Gaussian Elimination. I tried the following, rewriting the equations as:
$$ \left\{ \begin{array}{ll} \frac{\partial u (0,y)}{\partial x} = 0 = \frac{u_{1,j}-u_{-1,j}}{2\Delta x} \\ \frac{\partial u (x,10)}{\partial y} = 1 = \frac{u_{i,11}-u_{i,9}}{2\Delta y}\\ \end{array} \right. $$
Yielding ghost points:
$$ \left\{ \begin{array}{ll} u_{-1,j} = u_{1,j} \\ u_{i,11} = u_{i,9}+2 \Delta y\\ \end{array} \right. $$
$\Delta x \Delta x$ is defined as $h_2$. Subscript $j$ indicates the row ranging from $1$ to $M$, while subscript $i$ denotes the column, ranging from $1$ to $N$. I get the following equations for the $j^{th}$ row when writing down my equations when trying to solve the system via tridiagonal matrix: $$ \left\{ \begin{array}{ll} u_{-1,j} = u_{1,j} \\ \alpha_1u_{0,j} + \beta_1u{1,j} + \gamma_1u_{2,j}=f_{1,j}h_2 - u_{1,j-1} - u_{1,j+1}\\ \alpha_2u_{1,j} + \beta_2u{1,j} + \gamma_2u_{2,j}=f_{2,j}h_2 - u_{2,j-1} - u_{2,j+1}\\ \alpha_3u_{2,j} + \beta_3u{1,j} + \gamma_3u_{2,j}=f_{3,j}h_2 - u_{3,j-1} - u_{3,j+1}\\ \vdots\\ \alpha_{N-1}u_{N-2,j} + \beta_{N-1}u{1,j} + \gamma_{N-1}u_{2,j}=f_{N-1,j}h_2 - u_{N-1,j-1} - u_{N-1,j+1}\\ \alpha_{N}u_{N-1,j} + \beta_Nu{N,j} + \gamma_Nu_{N+1,j}=f_{N,j}h_2 - u_{N,j-1} - u_{N,j+1}\\ u_{i,11} = u_{i,9}+2 \Delta y\\ \end{array} \right. $$
I want to strictly solve the set of equations by using Gaussian Elimination, but I don't seem to find a way to eliminate $u_0$ in the 2nd equation of my set of equations. Is there a way to solve these equations strictly with this method?
It's not a tridiagonal matrix; most equations involve five unknowns. It is sparse, but having to flatten a matrix of unknowns into a vector of unknowns means we don't get the nice multidiagonal form that 1D Laplace equation has.
Also your system is incomplete: you have fewer equations than unknowns. You need to use the difference scheme at the points of Neumann boundary too, since these points have all four neighbors (some of them ghosts):
$$\alpha_1u_{-1,j} + \beta_1u_{0,j} + \gamma_1u_{1,j}=f_{0,j}h_2 - u_{0,j-1} - u_{0,j+1}$$
which I think means simply
$$ u_{-1,j} + u_{1,j} -4 u_{0,j} + u_{0,j-1} + u_{0,j+1} =0$$
Here $u_{-1,j}$ is $u_{1,j}$ by the boundary condition. Now you have a second equation with $u_{0,j}$.
That said, working Gaussian elimination by hand on this system is going to be very painful.