EigenFun[x_, y_, n_, m_] = 
      Cos[π*x/a*(n + 1/2)]*Cos[π*y/b*(m + 1/2)]

Coef[n_, m_] = 
      Integrate[EigenFun[x, y, n, m], {x, -a, a}, {y, -b, b}]/
        Integrate[(D[EigenFun[x, y, n, m], x, x] + 
                D[EigenFun[x, y, n, m], y, y])*
            EigenFun[x, y, n, m], {x, -a, a}, {y, -b, b}]

fii[x_, y_, N_, M_] := 
      Sum[Coef[i, j]*EigenFun[x, y, i, j], {i, 0, N}, {j, 0, M}]

Plot3D[
      fii[x, y, 5, 5] /. {a -> 1, b -> 1}, {x, -1, 
        1}, {y, -1, 1}]

LaplaceCoef[n_, 
        m_] = ((π*(n + 1/2)/a)^2 + (π*(m + 1/2)/
                    b)^2)*Coef[n, m]

LaplaceFii[x_, y_, N_, M_] := 
      Sum[LaplaceCoef[i, j]*EigenFun[x, y, i, j], {i, 0, N}, {j, 0, 
          M}]

Plot3D[
      LaplaceFii[x, y, 5, 5] /. {a -> 1, b -> 1}, {x, -1, 
        1}, {y, -1, 1}]

xlength = 1;
ylength = 1;
NIntegrate[((LaplaceFii[x, y, 2, 2] + 1) /. {a -> xlength, 
              b -> ylength})^2, {x, -xlength, 
        xlength}, {y, -ylength, ylength}]

xlength = 1;
ylength = 1;
Nmax = 10;
Mmax = 10;
L2norms = 
        Table[NIntegrate[((LaplaceFii[x, y, i, j] + 1) /. {a -> 
                      xlength, b -> ylength})^2, {x, -xlength, 
              xlength}, {y, -ylength, ylength}], {i, 0, Nmax}, {j, 0, 
            Mmax}];

ListPlot3D[L2norms]

DfiiDxCoef[n_, m_] = 
      Coef[n, m]*π/a*(n + 1/2)
DfiiDyCoef[n_, m_] = Coef[n, m]*π/b*(m + 1/2)

DfiiDxCoef[1, 1]

DfiiDx[x_, y_, N_, M_] := 
      Sum[-DfiiDxCoef[i, j]*Sin[π*x/a*(i + 1/2)]*
          Cos[π*y/b*(j + 1/2)], {i, 0, N}, {j, 0, 
          M}]
DfiiDy[x_, y_, N_, M_] := 
      Sum[-DfiiDyCoef[i, j]*Cos[π*x/a*(i + 1/2)]*
          Sin[π*y/b*(j + 1/2)], {i, 0, N}, {j, 0, M}]

Plot3D[
      DfiiDx[x, y, 5, 5] /. {a -> 1, b -> 1}, {x, -1, 
        1}, {y, -1, 1}]
Plot3D[
      DfiiDy[x, y, 5, 5] /. {a -> 1, b -> 1}, {x, -1, 
        1}, {y, -1, 1}]

apu[n_, m_] = 
      FullSimplify[
        Integrate[(-DfiiDxCoef[n, m]*Sin[π*x/a*(n + 1/2)]*
                  Cos[π*y/b*(m + 1/2)])^2 + (-DfiiDyCoef[n, 
                      m]*Cos[π*x/a*(n + 1/2)]*
                  Sin[π*y/b*(m + 1/2)])^2, {x, -a, 
            a}, {y, -b, b}]]

J[N_, M_] := 4*Sum[apu[i, j], {i, 0, N}, {j, 0, M}]

N[J[1, 1] /. {a -> 1, b -> 1}]

Nmax = 10;
Mmax = 10;
Jvalues = Table[J[i, j], {i, 0, Nmax}, {j, 0, Mmax}];

xlength = 1;
ylength = 1;
ListPlot3D[
      Jvalues /. {a -> xlength, b -> ylength}]

xlength = 20;
ylength = 1;
Nmax = 10;
Mmax = 10;
foppl = xlength*(ylength/2)^3/3;
ErrorValues = 
        Table[J[i, j] - foppl, {i, 0, Nmax}, {j, 0, 
            Mmax}];
ListPlot3D[
      Jvalues /. {a -> xlength, b -> ylength}]

xlength = 20;
ylength = 1;
Nmax = 10;
Mmax = 10;
foppl = 2*xlength*(ylength*2)^3/3;
ErrorValues = 
        Table[Abs[(J[i, j] - foppl)], {i, 0, Nmax}, {j, 0, 
            Mmax}];
ListPlot3D[
      ErrorValues /. {a -> xlength, b -> ylength}]