function ZP = potdiffGS(N,bc,M)
% POTDIFFGS  Solve potential eqn with finite diffs using 
%    Gauss-Seidel.  Compare the results and speeds from 
%    "potdiffGS(20,ones(1,19),M);" for M=2,5,10,20,40,80 to 
%    "potdiffS(20,ones(1,19));".
% (Ed Bueler; 2/13/02)
%
% uses finite differences to set up (N-1)^2 by (N-1)^2 matrix 
%    problem and solves it using the Gauss-Seidel method 
%    with M iterations

% see comments in potdiffS
G=numgrid('S',N+1); 
A=-delsq(G); 
b = zeros((N-1)^2,1);  
b(N-1:N-1:(N-1)^2)=-bc;

% Solve A x = b by Gauss-Seidel:  write as x = Tx + c and iterate
T=(1/4)*(A+4*eye((N-1)^2,(N-1)^2));
c=(-1/4)*b';
j=1;
Z=c; % start with c (saves a step over starting with zero)
while j<=M
   for k=1:(N-1)^2
      Z(k)= T(k,:)*Z'+c(k);
   end;
   j=j+1;
end;

% turn Z back into grid coordinates for display
ZP = zeros(N+1,N+1);
ZP(G>0)=Z(G(G>0)); % sneaky way of putting values into grid
ZP(N+1,:)=[0 bc 0];
ZP=ZP';

% plot
dx=1/N; xp=0:dx:1; yp=xp;
mesh(xp, yp, ZP);
xlabel('x axis'); ylabel('y axis'); axis([0 1 0 1 0 1.2]);
title('POTDIFFGS: Solution of potential eqn problem on unit square.');