Integrate over one variable of a 2D interpolating function returned from NDEigensystem

I’m trying to implement the answer for "Integrate only one variable of a 2D interpolating function" (https://mathematica.stackexchange.com/a/161962/73672) but for interpolating functions returned from NDEigenSystem it isn’t working.

First I solved for my required eigenfunctions and stored them as "funs":

ClearAll["Global`*"];
ClearAll[vals, funs, schröd];
A = 0.025;
Subscript[V, 0] = 1;
d = 2;
schröd = -A*d^2 π*D[ψ[n, φ], {φ, 2}] + 
   A/(4 π) (φ*φ*ψ[n, φ] + 
      2 I*D[ψ[n, φ], {n, 1}] - 
      D[ψ[n, φ], {n, 2}]) - 
   Subscript[V, 
     0] ((Cos[2 π*d*n]) + Cos[φ] - 20) ψ[
     n, φ];
Subscript[n, min] = -1/2; Subscript[n, max] = 1/2; 
Subscript[φ, min] = -π; 
Subscript[φ, max] = π;
Ω = 
  Rectangle[{Subscript[n, min], 
    Subscript[φ, min]}, {Subscript[n, max], 
    Subscript[φ, max]}];

{vals, funs} = 
  NDEigensystem[{schröd, 
    PeriodicBoundaryCondition[ψ[n, φ], 
     Subscript[φ, min] <= φ <= 
       Subscript[φ, max] && n == Subscript[n, max], 
     FindGeometricTransform[{{Subscript[n, min], 
         Subscript[φ, min]}, {Subscript[n, min], 
         Subscript[φ, max]}}, {{Subscript[n, max], 
         Subscript[φ, min]}, {Subscript[n, max], 
         Subscript[φ, max]}}][[2]]], 
    PeriodicBoundaryCondition[Exp[I 2 π n]*ψ[n, φ],
      Subscript[n, min] <= n <= Subscript[n, max] && φ == 
       Subscript[φ, max], 
     FindGeometricTransform[{{Subscript[n, min], 
         Subscript[φ, min]}, {Subscript[n, max], 
         Subscript[φ, min]}}, {{Subscript[n, min], 
         Subscript[φ, max]}, {Subscript[n, max], 
         Subscript[φ, 
          max]}}][[2]]]}, ψ, {n, φ} ∈ 
    Rectangle[{Subscript[n, min], 
      Subscript[φ, min]}, {Subscript[n, max], 
      Subscript[φ, max]}], 2];
Plot3D[{Evaluate[Abs[{funs[[1]][n, φ]  }^2]]}, {n, 
  Subscript[n, min], Subscript[n, max]}, {φ, 
  Subscript[φ, min], Subscript[φ, max]}, 
 PlotRange -> All, 
 PlotLabel -> 
  "Eig  fun    1 ", AxesLabel -> Automatic]

The answer to the other question was as follows:

da = 
  Flatten[
    Table[
      {t, tau, N @ Sin[2 (t + 3 tau)] Exp[-2 t - tau]}, 
      {t, 0, 2, 2/100}, {tau, 0, 5, 5/100}],
    1];
f = Interpolation @ da;
{{x1, x2}, {y1, y2}} = f["Domain"];
intx = Integrate[f[x, y], x] /. x -> x2;
    
nintx[y_?NumericQ] := Module[{x}, NIntegrate[f[x, y], {x, x1, x2}]];
    
Plot[nintx[y], {y, y1, y2}, PlotRange -> All]

However, trying to implement this myself gives an error after the third line, I believe because it reads n2 as 0.5:

f = funs[[1]];
{{n1, n2}, {φ1, φ2}} = f["Domain"];
intn = Integrate[f[n, φ], n] /. n -> n2;

General::ivar: 0.5 is not a valid variable.
Integrate::ilim: Invalid integration variable or limit(s) in 0.5.