You are almost there. All eigenvalues being zero means that the matrix is nilpotent: one of its powers is the zero matrix. In this case $H^3=0$. And since the error at $N$th step is controlled by $\sum_{n\ge N}\|H^n\|$, the error becomes zero at the third step.
Let's check this explicitly, writing $v=L_*^{-1}b$ for the vector added at each iteration:
$$\begin{align}x^{1}&=Hx^{0}+v \\
x^{2}&=Hx^{1}+C = H^2 x^0+Hv+v\\
x^{3}&=Hx^{2}+C = H^3 x^0+H^2v+Hv+v = \color{red}{H^2v+Hv+v} \\
x^{4}&=Hx^{3}+C = H^3v+H^2v+Hv+v = \color{red}{H^2v+Hv+v} \\
\dots
\end{align}$$
You may want check that the vector in red is indeed the solution, even thought the general theorem tells you so.
Ok, let's start with the method.
Gauss-Seidel iteration uses a decomposition of a matrix into a lower triangular matrix $L$ and a strictly upper triangular matrix $U$, $A = L + U$. Let's begin by defining our function and forming $A$, $L$, and $U$ based on a user's choice of $n$
function X = MyGaussSeidelExample(n)
A = 2*diag(ones(n,1)) - diag(ones(n-1,1),-1) - diag(ones(n-1,1),1); % We don't actually need to form this explicitly!
L = 2*diag(ones(n,1)) - diag(ones(n-1,1),-1);
U = -diag(ones(n-1,1),1)
b = 1/(n+1)^4*(1:n)'.^2; b(1) = b(1) + 1; b(n) = b(n) + 6; % Form b according to problem statement.
% (1:n) gives the vector [1 2 3 ... n]. (1:n)' makes it a column vector.
% (1:n)'.^2 squares each element. Pre-multiplying by 1/(n+1)^4 gives us the desired form.
% Manually add 1 to the first entry and 6 to the last.
x = zeros(n,1);
X = x;
end
This sets up your problem.
Now, the iteration rule is $Lx_{k+1} = b-Ux_k$. So, naively, we can solve $x_{k+1} = L^{-1}(b-Ux_k)$ at each step. Note, this won't end up being the best way to do it, but let's start with it.
function X = MyGaussSeidelExample(n)
A = 2*diag(ones(n,1)) - diag(ones(n-1,1),-1) - diag(ones(n-1,1),1); % We don't actually need to form this explicitly!
L = 2*diag(ones(n,1)) - diag(ones(n-1,1),-1);
U = -diag(ones(n-1,1),1);
b = 1/(n+1)^4*(1:n)'.^2; b(1) = b(1) + 1; b(n) = b(n) + 6; % Form b according to problem statement.
% (1:n) gives the vector [1 2 3 ... n]. (1:n)' makes it a column vector.
% (1:n)'.^2 squares each element. Pre-multiplying by 1/(n+1)^4 gives us the desired form.
% Manually add 1 to the first entry and 6 to the last.
x = zeros(n,1);
X = ones(n,1);
tol = 10e-8; % I assume you want an error tolerance to be small.
epsilon = 10*tol;
while (epsilon > tol)
x = inv(L)*(b-U*x);
epsilon = max(abs(X-x));
X = x;
end
total_error = max(abs(A\b - X)); % Compare the error to MATLAB's built in solver
end
Computing this in this manner using the inverse of $L$ is still very crappy. However, if you read the commentary here, it should be easy enough to directly compute each element of $x_{k+1}$ at each step without forming $L^{-1}$ directly. I'll leave that as an exercise for you. I have created the framework you need.
Best Answer
You can take the actual solution $x=x^{(k)}$ and get $$(L+D)x^{(k+1/2)} + Rx= b = (L+D+R)x.$$ If $L+D$ is invertible (it has to be, else the algorithm will not work), you get $x=x^{(k+1/2)}$. The same argument will now apply to the second equation, thus you will get $x=x^{(k+1)}$.
By subtracting the two equations, you can prove that if $x^{(k)}=x^{(k+1)}$ you have $Ax^{(k)}=b$.