Gauss-Seidel iteration in Matlab

Hi, I am trying to implement Gauss-Seidel in Matlab without using too many built in functions.

The way I am approaching it is to rewrite the system $\displaystyle Ax=b$ as

$\displaystyle x = D^{-1}(b-A_{off}\cdot x)$ where $\displaystyle A_{off}$ is the matrix $\displaystyle A$ with it's diagonal zeroed out.

Since I assume that the diagonal element of $\displaystyle A$ are non-zero, the elements in $\displaystyle D^{-1}$ are simply the reciprocals of the diagonal elements of $\displaystyle A$.

The first part of my code is just creating a strictly diagonally dominant matrix.

The loop converges to a solution, but it is not the right one, $\displaystyle Ax \neq b$. I have stared at it for quite some time now, and do not see what I have done wrong.

Code:

`format long;`

%-------------------------------TESTMATRIX------------------------------

n = 3;

A = zeros(n);

% The diagonal of A.

for i = 1:n

A(i,i) = 2*i;

end

% Superdiagonal

for i = 1:n-1

A(i,i+1)=i-1;

end

% Subdiagonal

for i = 2:n

A(i,i-1) = -i+1;

end

% Guessing a vector b in the system Ax=b.

b = zeros(1,n);

for i = 1:n

b(i) = i/2;

end

%-----------------------------TESTMATRIX END----------------------------

d = diag(A);

Aoff = A - diag(d);

x = zeros(1,n); % My first guess at a solution.

N = 10;

e = 10e-5;

while N >= 0

x_previous = x;

for i = 1:n

for j = 1:n

x(i) = x(i) + (Aoff(i,j)*x(j));

end

x(i) = (b(i)-x(i))/A(i,i);

end

if norm(x-x_previous,inf) < e

x

return

end

N = N - 1;

end