function ZP = potdiff(N)
%POTDIFF  Solve potential eqn with finite differences
% Try "potdiff(10);"
% (Ed Bueler; 1/18/02)
%
% Grid is N+1 by N+1 points in unit square (includes bdry).
% Interior N-1 by N-1 points have unknown values.
% Uses finite differences to set up (N-1)^2 by (N-1)^2 
%    matrix problem and solves it using Matlab's built-in "A\b" 
%    command (Gaussian elimination).
% There is dramatic increase in work with increase in N:
%    matrix A is (approx) N^2 by N^2
%    uses method (Gauss elim) which takes (N^2)^3=N^6 time.
% Grid indices are (j,k).  Global index is s.
% Compute global from grid by:  s = (k-1)(N-1)+j.

A = zeros((N-1)^2,(N-1)^2); b = zeros((N-1)^2,1);
dx=1/N; dy=dx;

for j=1:N-1
   for k=1:N-1
      s = (N-1)*(k-1)+j;
      sjp1 = s+1; % right neighbor
      sjm1 = s-1; % left neighbor
      skp1 = s+N-1; % upper neighbor
      skm1 = s-(N-1); % lower neighbor
      A(s,s) = -4;
      if j == 1 
         A(s,sjp1)=1;
         if k == 1
            A(s,skp1)=1;
         elseif k == N-1
            A(s,skm1)=1;
         else
            A(s,skp1)=1;
            A(s,skm1)=1;
         end;
      elseif j == N-1
         A(s,sjm1)=1;
         b(s)=-1;
         if k == 1
            A(s,skp1)=1;
         elseif k == N-1
            A(s,skm1)=1;
         else
            A(s,skp1)=1;
            A(s,skm1)=1;
         end;
      else
         A(s,sjp1)=1;
         A(s,sjm1)=1;
         if k == 1
            A(s,skp1)=1;
         elseif k == N-1
            A(s,skm1)=1;
         else
            A(s,skp1)=1;
            A(s,skm1)=1;
         end;
      end;
   end;
end;

% show A or use "spy"
% figure(2); spy(A); 
% title(['A is ' num2str((N-1)^2) ' by ' num2str((N-1)^2)]);
 
% solve A v = b --- faster program will use a sparse method
Z = A\b;

% turn Z back into grid coordinates
ZP = zeros(N+1,N+1);
for j=1:N-1
   for k=1:N-1
      s=(k-1)*(N-1)+j;
      ZP(j+1,k+1) = Z(s);
   end;
end;
ZP(N+1,:)=[0 ones(1,N-1) 0];
ZP=ZP';

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