function ZP = potfour(N,M)
%POTFOUR  Solve potential eqn with Fourier series
% Try "potfour(40,40);"
% (Ed Bueler; 1/18/02)
%
% Uses M terms in sum; work is proportional to N^2 * M.
% There is a maximum size for M of around 110, since  
%    sinh((2M-1)*pi) gets calculated, and it is a large number.

% create mesh of points
dx = 1/N;
xp=0:dx:1; yp=xp;
[XP,YP]=meshgrid(xp,yp);
%   This program is a good example of how "meshgrid" works.

% sum Fourier series solution
ZP=zeros(size(XP));
for j=M:-1:1  % count down for less rounding error
   JPI=(2*j-1)*pi;
   C=4/(JPI*sinh(JPI));
   ZP=ZP+C*sinh(JPI*XP).*sin(JPI*YP);
end;

% display
figure(1);
mesh(xp,yp,ZP);
xlabel('x axis'); ylabel('y axis');
axis([0 1 0 1 0 1.2]);
title('POTFOUR:  Soln of potential eqn problem on unit square.');