% 3D plot of solution to problem 3 on final, Math 421 Fall 2007 (Bueler)

NN=20; % number of terms in sum
x = linspace(0,1,201);    y = linspace(0,pi,31);
[xx,yy] = meshgrid(x,y);

zz = zeros(size(xx));
c = zeros(1,NN);
for k=0:NN-1
  n = 2 * k + 1;
  c(k+1) = -4 / (pi * n);
  % direct calculation easily overflows:
  %   ff = sinh(n * pi * (yy-pi)) ./ sinh(n * pi*pi));  
  small = exp(-2 * n * pi*pi);
  ff = ( exp(n*pi*(yy - 2 * pi)) - exp(-n*pi*yy) ) / (1 - small);
  zz = zz + c(k+1) * ff .* sin(n * pi * xx);
end

surf(x,y,zz,'LineStyle','none')
xlabel x, ylabel y, axis tight