Series[(1/4)*(F[x+h,y+h]+F[x-h,y+h]+F[x+h,y-h]+F[x-h,y-h]),{h,0,6}]

f[x_]:=x*(1-x);

f[x]

xmax=10;ymax=30;
const=1/10;ll=1;Clear[hh];

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*ll/xmax];hh[x,ymax]=0,{x,1,xmax-1}]

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

ll=1;
rel=10;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}]

nn=10;
Do[a[nnn]=2*NIntegrate[f[x]*Sin[nnn*x*Pi/ll],{x,0,ll}],{nnn,1,nn,2}];so=Sum[a[nnn]*Sin[nnn*Pi*x]*Exp[-nnn*Pi*y],{nnn,1,nn,2}]

xmax=10;ymax=30;
const=1/10;Clear[hh];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*ll/xmax];hh[x,ymax]=0,{x,1,xmax-1}];

rel=400;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,20]==0, bb=20;
Print[plo[qq]=Show[ListPlot[Table[{x,hh[x,bb]},{x,0,xmax}],PlotRange->All,PlotStyle->PointSize[0.02],DisplayFunction->Identity],Plot[(so/.{y->(ll/xmax)*bb,x->(ll/xmax)*xx}),{xx,0,xmax},DisplayFunction->Identity],
DisplayFunction->$DisplayFunction]];
(*ListPlot3D[Table[hh[x,y],{x,0,xmax},{y,0,ymax}],PlotRange->All]*)],{qq,1,rel}]

Show[Table[plo[qq],{qq,10,100,10}],spl]

spl=Plot[(so/.{y->(ll/xmax)*bb,x->(ll/xmax)*xx}),{xx,0,xmax}]

const=1/10;Clear[hh];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*ll/xmax];hh[x,ymax]=0,{x,1,xmax-1}];rel=400;Do[
Do[hhnew[x,y]=(4/20)*(hh[x+1,y]+hh[x-1,y]+hh[x,y+1]+hh[x,y-1])+(1/20)*(hh[x+1,y+1]+hh[x-1,y+1]+hh[x+1,y-1]+hh[x-1,y-1])
,{x,1,xmax-1},{y,1,ymax-1}];
Do[hh[x,y]=hhnew[x,y],{x,1,xmax-1},{y,1,ymax-1}];bb=20;
If[Mod[qq,50]==0,Show[ListPlot[Table[{x,hh[x,bb]},{x,0,xmax}],PlotRange->All,PlotStyle->PointSize[0.02],DisplayFunction->Identity],Plot[(so/.{y->(ll/xmax)*bb,x->(ll/xmax)*xx}),{xx,0,xmax},DisplayFunction->Identity],
DisplayFunction->$DisplayFunction]],{qq,1,rel}]