function gs5ptdemo
%GS5PTDEMO  graphics demonstration of convergence of Gauss-Seidel
%
% demonstration of Gauss-Seidel for 5-point Laplacian, illustrating
% the smoothing property
% Douglas N. Arnold, 2009-10-10

niters = 250;

% supply f
%f = @(x,y) -2*pi^2*sin(pi*x).*sin(pi*y);
f = @(x,y) 190*x.^2.*(1-y)-26.4*exp((x-.25).^2+(y-.25).^2);

% set mesh size
n = 64;
h = 1/n;
m = n-1;
x = linspace(0,1,n+1);   % x grid points including boundaries
y = linspace(0,1,n+1);   % y grid points including boundaries

[X,Y] = meshgrid(x,y);      % 2d arrays of x,y values
X = X';                     % transpose so that X(i,j),Y(i,j) are
Y = Y';                     % coordinates of (i,j) grid point, i.e., ((i-1)h,(j-1)h)

Iint = 2:n;                 % indices of interior points in x
Jint = 2:n;                 % indices of interior points in y
Xint = X(Iint,Jint);        % interior points
Yint = Y(Iint,Jint);

fh = f(Xint,Yint);          % fh is a mesh function on interior grid points

%%%%%%  Direct solve %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% set uh on the boundary to given values (interior values will be ignored)
uh = zeros(n+1,n+1);        % can use true solution if it is known

% form 5-point Laplacian matrix (times h^2)
I = speye(m);
e = ones(m,1);
A = spdiags([e,-4*e,e],[-1,0,1],m,m);
J = spdiags([e,e],[-1,1],m,m);
Lh = kron(I,A) + kron(J,I);

fh = f(Xint,Yint);          % fh is a mesh function on interior grid points
% multiply right hand side by h^2 and adjust for boundary condition
gh = h^2*fh;
gh(:,1) = gh(:,1) - uh(Iint,1);
gh(:,m) = gh(:,m) - uh(Iint,n+1);
gh(1,:) = gh(1,:) - uh(1,Jint);
gh(m,:) = gh(m,:) - uh(n+1,Jint);

% convert the right hand side from a grid function into a vector
F = reshape(gh,m*m,1);

% Solve the linear system for the grid function uh on the interior, written
% as vector
U = Lh\F;  

% convert from vector into grid function use this for uh at the interior
% grid points (uh is already given at the boundary grid points)
uh(Iint,Jint) = reshape(U,m,m);
% plot results:
figure(1)
clf
uhplot=surf(X,Y,uh);
set(uhplot,'facecolor','yellow')
set(gca,'zlim',[-1,1])
l=light('position',[-5,5,10]);
lighting gouraud
%%%%%%  End direct solve %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%% Gauss-Seidel iteration %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% set up initial iterate (random grid function vanishing on boundary)
u = zeros(n+1,n+1);
rand('state',0); u(Iint,Jint) = 2*rand(m,m)-1;

figure(2)
clf
uplot=surf(X,Y,u);
set(uplot,'facecolor','yellow')
set(gca,'zlim',[-1,1])
hold on
l=light('position',[-5,5,10]);
lighting gouraud
% print max norm error
fprintf('iteration    0.   max norm error: %9.5f\n',max(max(abs(uh-u))))

% iterate, updating plot each time
for niter = 1:niters
  pause
  for i=2:n
    for j=2:n
      u(i,j) = (u(i+1,j)+u(i-1,j)+u(i,j+1)+u(i,j-1)-h^2*fh(i-1,j-1))/4;
    end
  end
  clf
  uplot=surf(X,Y,u);
  set(uplot,'facecolor','yellow')
  set(gca,'zlim',[-1,1])
  l=light('position',[-5,5,10]);
  lighting gouraud
  % print max norm error
  fprintf('iteration %4d.   max norm error: %9.5f\n',niter,max(max(abs(uh-u))))
end