% Example 5.3a (Physics 330)

clear;close all;

% set some constants

a=2;b=1;Vo=1;

% build the x and y grids

Nx=80;Ny=40;
dx=2*b/Nx;dy=a/Ny;
x=-b:dx:b;y=0:dy:a;

% build the 2-d grids for plotting

[X,Y]=meshgrid(x,y);

% set the number of terms to keep
% and do the sum

Nterms=20;

% zero V out so we can add into it

V=zeros(Ny+1,Nx+1);

% add the terms of the sum into V

for m=0:Nterms
   V=V+cosh((2*m+1)*pi*X/a)/cosh((2*m+1)*pi*b/a).*sin((2*m+1)*pi*Y/a)/(2*m+1);
end

% put on the multiplicative constant

V=4*Vo/pi*V;

% surface plot the result

surf(X,Y,V)
xlabel('x');
ylabel('y');
zlabel('V(x,y)')