clear
N=64;h=1/(N+1); %domain [0,1]*[0,1]
xx=h:h:1-h;[X,Y] = ndgrid(xx,xx); %interior uniform mesh
ffun=@(x,y) -2*pi^2*(cos(2*pi*x).*sin(pi*y).^2+sin(pi*x).^2.*cos(2*pi*y));
F=h^2*ffun(X,Y); b=F(:); %right-hand-side function values as a column vector
A=gallery('poisson',N); %generate SPD matrix A in (3.19)
u=A\b; %solve by sparse direct solver, most accurate
%U=reshape(u,N,N); mesh(X,Y,U) %reshape and plot the 2D solution u(x,y)

omg=(1.5:0.01:2);
tol=1e-7; itlist=[]; errlist=[];
for w=omg
[x,iter]=SOR(A,b,w,tol);
itlist=[itlist;iter];
errlist=[errlist;norm(A*x-b)/norm(b)];
end
subplot(2,1,1)
plot(omg,itlist); hold on
omg_opt=2/(1+sin(pi*h))
[minit,idx]=min(itlist);
omg_app=omg(idx)
plot(omg_opt,itlist(idx),'*')
xlabel('{\omega}'); ylabel('Iteration numbers');
subplot(2,1,2)
semilogy(omg,errlist)
xlabel('{\omega}'); ylabel('Relative residual norm');

function [x,iter]=SOR(A,b,omega,tol)
%Input: A is the system matrix
% b is the right hand side vector
% omega is the relaxation factor used in SOR
% tol is the stopping tolerance for relative residual reduction
%Output: x is solution of Ax=b, iter is the used iteration numbers

D = diag(diag(A)); % diagonal part of A as matrix
L = tril(A); % lower triangular part of A
Eomega = ((1-omega)*D + omega*L)/omega; % for SOR

x = zeros(size(b)); %initial guess
r = b - A*x; r0=norm(r);
for iter=1:1000
x = x + Eomega \ r;
r = b - A*x;
if norm(r)/r0 < tol %relative residual
break;
end
end
end