%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Demonstration program for solution of 
% 
%  u_xx = f(x)
%
% using O(h^2) finite difference tecnhiques and successive overrelaxation
% (or Gauss-Seidel relaxation for omega = 1.0).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
level = 5;
nx = 2^level + 1;
x = linspace(0, 1, nx);
h = x(2) - x(1);
f = -pi^2 * sin(pi * x);
u = zeros(1, nx);
% omega = 1      -> Gauss-Seidel 
% 1 < omega <= 2 -> SOR (Successive Overrelaxation)
omega = 1.8;
clf;
hold on;

epsilon = 1.0e-4;
du = 1;
sweep = 0;
while du > epsilon
   sweep = sweep + 1;
   du = 0;
   for j = 2 : nx - 1
      uold = u(j);
      ugs = 0.5 * (u(j+1) + u(j-1) - h^2*f(j));
      u(j) = omega * ugs + (1 - omega) * uold;
      du = du + (uold - u(j))^2;
   end
   du = sqrt(du/nx);
   if rem(sweep,50) == 0
      fprintf('Norm of solution change after %d sweeps is %g\n', ...
         sweep, du);
      plot(x, u); plot(x,u,'or');
   end
end
