Diffusion define a function on 0 < x < 2, 0 < y < 3 restart; with(CodeGeneration): phi := (x,y) -> (cos(Pi/2 * x))^5 * (sin(Pi/3 * y))^3 + cos(Pi/2 * x)^2*sin(Pi/3 * y)/2; plot plot3d(phi,0..2,0..3,grid=[40,40],axes=boxed,labels=[x-axis,y-axis,z-axis],numpoints=300); calculate fourier coefficients, a(k,l) a:= (k,l) -> value(integrate(integrate(phi(x,y)*cos(k*Pi*x/2)*sin(l*Pi*y/3),x=0..2),y=0..3)); calculate normalization, nz(k,l) nz:= (k,l) -> value(integrate(integrate((cos(k*Pi*x/2)*sin(l*Pi*y/3))^2,x=0..2),y=0..3)); for k from 0 to 5 do for l from 1 to 3 do b[k,l] := a(k,l)/nz(k,l) end do end do; for l from 1 to 3 do seq(b[k,l],k=0..5) end do; check parseval blist2 := [seq(seq((b[k,l])^2*(nz(k,l)/(3/2)),k=0..5),l=1..3)]: sum(blist2[jj],jj=1..18)*(2/2)*(3/2); value(int(int((phi(x,y))^2,x=0..2),y=0..3)); plot fourier representation flist := [seq(seq(b[k,l]*cos(k*Pi*x/2)*sin(l*Pi*y/3),k=0..5),l=1..3)]: Phi := (x,y) -> value(sum(flist[jj],jj=1..18)); plot3d(Phi(x,y),x=0..2,y=0..3,grid=[40,40],axes=boxed,labels=[x-axis,y-axis,z-axis],numpoints=300); plot error from original function plot3d(Phi(x,y)-phi(x,y),x=0..2,y=0..3,grid=[40,40],axes=boxed,labels=[x-axis,y-axis,z-axis],title=['Phi(x,y)-phi(x,y)'],numpoints=300); calculate decay coefficients, dk(k,l) for k from 0 to 5 do for l from 1 to 3 do dk[k,l] := (k*Pi/2)^2 + (l*Pi/3)^2 end do end do; dklist := [seq(seq(exp(-dk[k,l]*t),k=0..5),l=1..3)]: animation plot, u(x,y,t) u := (x,y,t) -> value(sum(dklist[jj]*flist[jj],jj=1..18)); plots[animate](plot3d,[u(x,y,t),x=0..2,y=0..3],t=0..1,frames=100);