differential equations - How to prevent instability blow up in NDSolve?

I have the following code to solve a PDE:

e = 2.5;
xmax = 5;
ymax = 5;
sol[x_, y_] = f[x, y] /. First@NDSolve[{
     -D[f[x, y], x, x] - D[f[x, y], y, y] == e f[x, y],
     Derivative[0, 1][f][x, -ymax] == Cos[\[Pi]/(2 xmax) x],
     f[x, -ymax] == 0,

     f[-xmax, y] == 0,
     f[xmax, y] == 0
     }, f[x, y], {x, -xmax, xmax}, {y, -ymax, ymax}]

I'd expect to get a Sin[]-like solution in y direction, but what I get instead is this:

Plot3D[sol[x, y], {x, -xmax, xmax}, {y, -ymax, ymax}, PlotRange -> {-1, 1}, AxesLabel -> {"x", "y", "f"},  MaxRecursion -> 3]

If I try using something like MaxStepSize -> 0.25, then the blow-up is just more frequent in x direction, and becomes visible at even smaller y values:

What can I do to prevent this blow-up and make NDSolve give the expected solution?

EDIT

In fact, the problem above is a minimal example showing the instability. What I'd actually like to solve is the equation with additional +U[x,y]f[x,y] on the LHS, where U[x,y]=-70Exp[-x^2-y^2]. I.e. the equation would look like

-D[f[x, y],x,x]-D[f[x,y],y,y]-70Exp[-x^2-y^2]f[x,y]==e f[x,y]

So, the answer I'm looking for should be extensible to this case. The answer by @xzczd solves the problem with minimal example, but unfortunately fails to extend to this one.