how can I change the Gauss-Seidel method to SOR methodcode in Matlab?

The question has shows that In implementing SOR methodin MATLAB, one should not calculate Tw and cw by formulas Tw = (D-wL)^(-1)[(1-w)D+wU)] and Cw = w(D-wL)^(-1)b , where w stands foromega and using MATLAB's built-in inv function, since this functionrequires O(n^3) flops and therefore the
whole matter loses its point.

I have tried for many times but I can't get the correctanswers. Please help Thanks :)

Execute code:

function x = GaussSeidel(A,b,x0,tol,kmax,output)
% This function returns an approximate solution of a
% linear system A*x=b, obtained by the Gauss-Seidel method.
%
% x0 is an initial approximation,
% tol is tolerance,
% kmax is the maximum number of iterations,
% output is a parameter which regulates displays:
% 0 - no display,
% 1 - some display,
% 2 - detailed display.
%
% Input: A, n by n matrix,
% b, n by 1 column,
% x0, initial approximation,
% tol, tolerance,
% kmax, maximum number of iterations,
% output, a parameter which regulates displays.
%
% Output: x, a solution.

[m,n]=size(A);
if m~=n,
error('A is not square');
end
m=length(b);
if m~=n,
error('Dimensions of A and B do not agree');
end
% Write A as D-L-U
[D,L,U] = DLU_decomposition(A);
% Calculate the matrix T and the vector c
T = L+U;
c=b;
for i=1:n,
T(i,:) = T(i,:)/D(i,i); % T = D^(-1)*(L+U)
c(i) = b(i)/D(i,i); % c = D^(-1)*b
end
% Calculate x
x=x0; % initial approximation
if output >= 2
disp(' k TOL x1 x2');
end
for k=1:kmax,
xprev=x;
for i=1:n
x(i)=c(i);
for j=1:n
if j~=i
x(i)=x(i)+T(i,j)*x(j);
end
end
end
TOL=norm(x-xprev,inf);
if output >= 2
out=[k, TOL, x(1), x(2)];
disp(out);
end
if TOLif output >= 1
disp('The Gauss-Seidel method has converged.');
end
break
end
end
if output >= 1
if TOL>=tol
disp('The Gauss-Seidel method did not converge.');
end
s=sprintf('The norm of residual vector A*x-b is%e.',norm(A*x-b));
disp(s);
end
end