% setup n and h
n=20;
h=1/n;

% setup A.
A=diag(-2*ones(n-1,1))+diag(ones(n-2,1),-1)+diag(ones(n-2,1),1);
A=A/h^2;

%setup b. Here, we try to solve u''=-sin(Pi*x), u(0)=u(1)=0.
x=(h:h:1-h/2)';
b=-sin(pi*x);

%the exact solution xexact
uexact=A\b;

%Initial guess for Ax=b
rand('state',0);
u=(rand(n-1,1)-0.5*ones(n-1,1));

%Now solve the problem using GS method.
for i=1:400,
    plot(x,uexact-u);
    pause(0.1);
    unew=SOR(2/3,n,h,u,b);
    u=unew;
end