This maple worksheet produce Fourier sine expansions, plots, and approximations of functions of two variables.
> with(plots);

[animate, animate3d, changecoords, complexplot, complexplot3d,

    conformal, contourplot, contourplot3d, coordplot, coordplot3d,

    cylinderplot, densityplot, display, display3d, fieldplot,

    fieldplot3d, gradplot, gradplot3d, implicitplot, implicitplot3d,

    inequal, listcontplot, listcontplot3d, listdensityplot, listplot,

    listplot3d, loglogplot, logplot, matrixplot, odeplot, pareto,

    pointplot, pointplot3d, polarplot, polygonplot, polygonplot3d,

    polyhedraplot, replot, rootlocus, semilogplot, setoptions,

    setoptions3d, spacecurve, sparsematrixplot, sphereplot, surfdata,

    textplot, textplot3d, tubeplot]

visualize the mode corresponding to m=K and n=L
> K:=2;
> L:=2;

> plot3d(sin(K*x)*sin(L*y),x=0..Pi,y=0..Pi);
> animate3d(sin(sqrt(K^2+L^2)*t)*sin(K*x)*sin(L*y),x=0..Pi,y=0..Pi,t=0..2*Pi,frames=20);

                                K := 2


                                L := 2

input the numbers N  and M of terms  to be used in the  partial sums 
> M:=5;

                                M := 5

> N:=3;

                                N := 3

input the initial position  f(x,y)  and initial velocity g(x,y)
> f:=(x,y)->(x-Pi)*sin(x/2)*sin(y);
> g:=(x,y)->0;

              f := (x, y) -> (x - Pi) sin(1/2 x) sin(y)


                                g := 0

define two  matrix of coefficients 
> a:=array(1..M,1..N);
> b:=array(1..M,1..N);

                    a := array(1 .. 5, 1 .. 3, [])


                    b := array(1 .. 5, 1 .. 3, [])

compute the Fourier coefficients of the function f(x,y) and print the result as a matrix
> for n to N do for m to M do a[m,n]:=2/Pi*int(sin(n*y)*2/Pi*int(sin(m*x)*f(x,y),x=0..Pi),y=0..Pi) od: od:
> print(a);

                       [        1            ]
                       [- 32/9 ----    0    0]
                       [        Pi           ]
                       [                     ]
                       [  64   1             ]
                       [- --- ----     0    0]
                       [  225  Pi            ]
                       [                     ]
                       [   96   1            ]
                       [- ---- ----    0    0]
                       [  1225  Pi           ]
                       [                     ]
                       [  128   1            ]
                       [- ---- ----    0    0]
                       [  3969  Pi           ]
                       [                     ]
                       [  160   1            ]
                       [- ---- ----    0    0]
                       [  9801  Pi           ]

compute the Fourier coefficients of the function g(x,y) and print the result as a matrix

> for n to N do for m to M do b[m,n]:=2/Pi*int(sin(n*y)*2/Pi*int(sin(m*x)*g(x,y),x=0..Pi),y=0..Pi) od: od:
> print(b);



                            [0    0    0]
                            [           ]
                            [0    0    0]
                            [           ]
                            [0    0    0]
                            [           ]
                            [0    0    0]
                            [           ]
                            [0    0    0]

write a double Fourier series for f(x,y) and g(x,y)
> sf:=(x,y)->sum('sum('a[j,k]*sin(j*x)','j'=1..M)*sin(k*y)','k'=1..N);

> sg:=(x,y)->sum('sum('b[j,k]*sin(j*x)','j'=1..M)*sin(k*y)','k'=1..N);

                      N     /   M                      \
                    -----   | -----                    |
                     \      |  \                       |
   sf := (x, y) ->    )    '|   )    'a[j, k] sin(j x)'| sin(k y)'
                     /      |  /                       |
                    -----   | -----                    |
                   'k' = 1  \'j' = 1                   /


                      N     /   M                      \
                    -----   | -----                    |
                     \      |  \                       |
   sg := (x, y) ->    )    '|   )    'b[j, k] sin(j x)'| sin(k y)'
                     /      |  /                       |
                    -----   | -----                    |
                   'k' = 1  \'j' = 1                   /

plot the function f(x,y) and its Fourier series approximation
> plot3d(f(x,y),x=0..Pi,y=0..Pi,title=`f(x,y)`);
> plot3d(sf(x,y),x=0..Pi,y=0..Pi,title=`Partial sum approximation of f(x,y)`);
plot the function g(x,y) and its Fourier series approximation
> plot3d(g(x,y),x=0..Pi,y=0..Pi,title=`g(x,y)`);
> plot3d(sg(x,y),x=0..Pi,y=0..Pi,title=`Partial sum approximation of g(x,y)`);
write the (partial) Fourier series solution
> u:=(x,y,t)->sum('sum('(a[m,n]*cos(sqrt(m^2+n^2)*t)+b[m,n]*(n^2+m^2)^(-1/2)*sin(sqrt(m^2+n^2)*t))*sin(m*x)','m'=1..M)*sin(n*y)','n'=1..N);


                     N     /   M
                   -----   | -----
                    \      |  \
u := (x, y, t) ->    )    '|   )    '
                    /      |  /
                   -----   | -----
                  'n' = 1  \'m' = 1

    /                                                 2    2    \
    |                  2    2       b[m, n] sin(sqrt(m  + n ) t)|
    |a[m, n] cos(sqrt(m  + n ) t) + ----------------------------|
    |                                         2    2 1/2        |
    \                                       (m  + n )           /

             \
             |
             |
    sin(m x)'| sin(n y)'
             |
             |
             /

plot the solution for a particular value of t (for example t=Pi) and animate the solution
> plot3d(u(x,y,Pi),x=0..Pi, y=0..Pi,title=`Solution at t=Pi`);
> animate3d(u(x,y,t),x=0..Pi, y=0..Pi, t=0..2*Pi, frames=20);
the solution to a problem we did in class
> p:=(x,y,t)->1/2*t-1/4*sin(2*t)*cos(2*x);
> animate3d(p(x,y,t),x=0..Pi, y=0..Pi, t=0..2*Pi, frames=20);

           p := (x, y, t) -> 1/2 t - 1/4 sin(2 t) cos(2 x)

>