(* this is primitive program to check 2.10 in Jackson; *)

f[x_]:=1;
xmax=10;ymax=30;
const=1/10;Clear[hh,ll,so];
nn=10;
Do[a[nnn]=(2/xmax)*NIntegrate[f[x]*Sin[nnn*x*Pi/xmax],{x,0,xmax}],{nnn,1,nn,2}];so[x_,y_]:=Sum[a[nnn]*Sin[nnn*Pi*x/xmax]*Exp[-nnn*Pi*y/xmax],{nnn,1,nn,2}];so[x,y]

ListPlot3D[Table[so[x,y],{x,0,xmax},{y,0,ymax}],PlotRange->All]

Do[hh[x,y]=const,{x,1,xmax-1},{y,1,ymax-1}];Do[hh[0,y]=0;hh[xmax,y]=0,{y,0,ymax}];Do[hh[x,0]=f[x/xmax];hh[x,ymax]=0,{x,1,xmax-1}]

ListPlot3D[Table[hh[x,y],{x,0,xmax},{y,0,ymax}],PlotRange->All]

Do[hh[x,y]=const,{x,1,xmax-1},{y,1,ymax-1}];Do[hh[0,y]=0;hh[xmax,y]=0,{y,0,ymax}];Do[hh[x,0]=f[x/xmax];hh[x,ymax]=0,{x,1,xmax-1}]

rel=3;Do[
Do[hh[2*x,2*y]=(1/4)*(hh[2*x+1,2*y]+hh[2*x-1,2*y]+hh[2*x,2*y+1]+hh[2*x,2*y-1]),{x,1,xmax/2-1},{y,1,ymax/2-1}];

Do[hh[2*x-1,2*y-1]=(1/4)*(hh[2*x-1,2*y]+hh[2*x-1,2*y-2]+hh[2*x,2*y-1]+hh[2*x-2,2*y-1]),{x,1,xmax/2},{y,1,ymax/2}];Do[hh[2*x,2*y-1]=(1/4)*(hh[2*x-1,2*y-1]+hh[2*x+1,2*y-1]+hh[2*x,2*y]+hh[2*x,2*y-2]),{x,1,xmax/2-1},{y,1,ymax/2}];Do[hh[2*x-1,2*y]=(1/4)*(hh[2*x-1,2*y+1]+hh[2*x-1,2*y-1]+hh[2*x,2*y]+hh[2*x-2,2*y]),{x,1,xmax/2},{y,1,ymax/2-1}];
If[Mod[qq,1]==0,Print[ListPlot3D[Table[hh[x,y],{x,0,xmax},{y,0,ymax}],PlotRange->All]]],{qq,1,rel}]

yfixed=1
Do[hh[x,y]=const,{x,1,xmax-1},{y,1,ymax-1}];Do[hh[0,y]=0;hh[xmax,y]=0,{y,0,ymax}];Do[hh[x,0]=f[x/xmax];hh[x,ymax]=0,{x,1,xmax-1}]

rel=100;Do[
(*even-even*)
Do[hh[2*x,2*y]=(1/4)*(hh[2*x+1,2*y]+hh[2*x-1,2*y]+hh[2*x,2*y+1]+hh[2*x,2*y-1]),{x,1,xmax/2-1},{y,1,ymax/2-1}];
(*odd-odd*)
Do[hh[2*x-1,2*y-1]=(1/4)*(hh[2*x-1,2*y]+hh[2*x-1,2*y-2]+hh[2*x,2*y-1]+hh[2*x-2,2*y-1]),{x,1,xmax/2},{y,1,ymax/2}];Do[hh[2*x,2*y-1]=(1/4)*(hh[2*x-1,2*y-1]+hh[2*x+1,2*y-1]+hh[2*x,2*y]+hh[2*x,2*y-2]),{x,1,xmax/2-1},{y,1,ymax/2}];Do[hh[2*x-1,2*y]=(1/4)*(hh[2*x-1,2*y+1]+hh[2*x-1,2*y-1]+hh[2*x,2*y]+hh[2*x-2,2*y]),{x,1,xmax/2},{y,1,ymax/2-1}];
sec[qq]=Table[{x,hh[x,yfixed]},{x,0,xmax}];
If[Mod[qq,30]==0,Print[ListPlot3D[Table[hh[x,y],{x,0,xmax},{y,0,ymax}],PlotRange->All]]],{qq,1,rel}]

relax=ListPlot[Table[{x,hh[x,yfixed]},{x,0,xmax}],PlotStyle->{RGBColor[1,0,0],PointSize[0.02]}]

fourier=ListPlot[Table[{x,so[x,yfixed]},{x,0,xmax}],PlotStyle->{RGBColor[0,0,1],PointSize[0.02]}]

exact=Plot[(2/Pi)*ArcTan[Sin[Pi*x/xmax]/Sinh[Pi*yfixed/xmax]],{x,0,xmax}]

Show[exact,fourier]

Show[exact,relax]