nmax = 50;
L = 1.;

y = double(0.:0.01:1.);
x = double(0.:0.01:1.);

Pi = double(pi)


for i = 1:length(x);
  for j = 1:length(y);
  u1(i,j) = 0.;
for n = 1:2:nmax;
term = sin(n*Pi*x(i)/L) * (-1)^((n-1)/2) / (n*Pi/L);
u1(i,j) = u1(i,j) - term * ( sinh(n*Pi/2.) * sinh(n*Pi*y(j)/L) / sinh(n*Pi));
if (y(j) >= L/2.) 
  u1(i,j) = u1(i,j) + term* sinh(n*Pi*(y(j)-L/2.)/L) ;
end
end
end
end

for i = 1:length(x);
  for j = 1:length(y);
  u2(i,j) = 0.;
for n = 1:2:nmax;
term = sin(n*Pi*y(j)/L) * (-1)^((n-1)/2) / (n*Pi/L);
u2(i,j) = u2(i,j) - term * ( sinh(n*Pi/2.) * sinh(n*Pi*x(i)/L) / sinh(n*Pi));
if (x(i) >= L/2.) 
  u2(i,j) = u2(i,j) + term* sinh(n*Pi*(x(i)-L/2.)/L) ;
end
end
end
end

surf(x,y,-(u1+u2)/2);



