% Solution of -u'' = f, u(0)=u(1)=0 by matrix iteration
% 
%     Let v = approximate solution 

N=20;
h=1/N;
x=linspace(0,1,N+1)';
x1=x(2:N);
D = 2*diag(diag(ones(N-1)))/h^2;
U= diag(diag(ones(N-2)),-1)/h^2;
L=diag(diag(ones(N-2)),1)/h^2;
A=D-U-L;
v0=sin(x1*pi);

% Jacobi Iteration
figure(1)
plot(x1,v0)
title('Jacobi Iteration')
xlabel('x')
ylabel('u')
hold on
vold=v0;
P = D^(-1)*(L+U);
for i=1:100
   v=P*vold;
   vold=v;
   plot(x1,vold)
   errJ(i) = max(v);
end
hold off
figure(2)
subplot(2,1,2)
plot(errJ)
title('Jacobi Iteration')
xlabel('Iteration')
ylabel('Error')

% Weighted Jacobi Iteration
figure(3)
plot(x1,v0)
title('Weighted Jacobi Iteration')
xlabel('x')
ylabel('u')
hold on
vold=v0;
P = D^(-1)*(L+U);
w = 2/3;
PWJ=(w*P+(1-w)*eye(N-1));
for i=1:100
   v=PWJ*vold;
   vold=v;
   plot(x1,vold)
   errWJ(i) = max(v);
end
hold off
figure(2)
subplot(2,1,1)
plot(errWJ)
title('Weighted Jacobi Iteration')
xlabel('Iteration')
ylabel('Error')

% Gauss-Seidel Iteration
figure(4)
plot(x1,v0)
title('Gauss-Seidel Iteration')
xlabel('x')
ylabel('u')
hold on
vold=v0;
P = (D-L)^(-1)*(U);
for i=1:100
   v=P*vold;
   vold=v;
   plot(x1,vold)
   errG(i) = max(v);
end
hold off
figure(2)
subplot(2,1,1)
plot(errG)
title('Gauss-Seidel Iteration')
xlabel('Iteration')
ylabel('Error')

