restart:

with(plots):

Vibrating square membrane

Initial conditions

psi0:=0;

v0:=2*exp(-(x^2+y^2)*4^2);

Normal frequencies

for m from 1 to 20 do

for n from 1 to 20 do

lambda[m,n] := Pi*sqrt(n^2+m^2);

od;

Normal mode coefficients for displacement

for m from 1 to 20 do

for n from 1 to 20 do

A[m,n] := evalf(Int(Int(psi0*sin(n*Pi*x)*sin(m*Pi*y),x=0..1),y=0..1)*4);

od;

Normal mode coefficients for velocity

for m from 1 to 20 do

for n from 1 to 20 do

B[m,n] := evalf(Int(Int(v0*sin(n*Pi*x)*sin(m*Pi*y),x=0..1),y=0..1)*4/lambda[m,n]);

od;

Solution

psi := add(add(sin(n*Pi*x)*sin(m*Pi*y)*(A[m,n]*cos(lambda[m,n]*t) + B[m,n]*sin(lambda[m,n]*t)),m=1..10),n=1..10):

animate3d([x,y,psi], x=0..1, y=0..1, t=0..20, frames=21, scaling=constrained, orientation=[45,65],lightmodel=light1, style=patchnogrid);

Nodal lines

nodes := (n,m)->sin(n*Pi*x)*sin(m*Pi*y) + sin(m*Pi*x)*sin(n*Pi*y);

contourplot(nodes(4,8),x=0..1,y=0..1,contours=[0],grid=[200,200]);