Here is the Mathematica code modeling of 2D vibrations of the rectangular region excited by a sine function.

k[x_, y_] = Sin[x + y]

Sin[x + y]

a[n_, m_] =

4 / (l1 l2) Integrate[

Integrate[

k[x, y] Sin[n Pi x / l1] Sin[m Pi y / l2], {x, 0, l1}], {y, 0, l2}]

-(4 (-2 n \[Pi] Cos[2] Sin[m \[Pi]] +

n \[Pi] Cos[n \[Pi]] (m \[Pi] Sin[4] + 2 Cos[6] Sin[m \[Pi]]) –

4 m \[Pi] Cos[4] Sin[n \[Pi]] +

8 Sin[6] Sin[m \[Pi]] Sin[n \[Pi]] +

m \[Pi] Cos[

m \[Pi]] (n \[Pi] Sin[2] – n \[Pi] Cos[n \[Pi]] Sin[6] +

4 Cos[6] Sin[n \[Pi]])))/((-4 + m^2 \[Pi]^2) (-16 +

n^2 \[Pi]^2))

Manipulate[

Plot3D[Sum[

Sum[a[n, m] Cos[

c Sqrt[((n Pi x / l1)^2) + ((m Pi y / l2)^2)] t] Sin[

n Pi x / l1] Sin[m Pi y / l2], {n, 1, B}], {m, 1, M}], {x, 0,

10}, {y, 0, 15}, PlotRange -> {-10, 10}]

,

{t, 0, 10},

{B, 1, 10},

{M, 1, 10},

{l1, Exp[-25], 12},

{l2, Exp[-25], 10},

{c, 0, 5}

]