qx = 0; b = 2; ycutoff = 15.0; Nless = 50;
lhs = (b^2 (-2 + qx) y (3 + (-2 + qx) y) - 2 b (-2 + qx) (1 + 2 qx y) + 4 qx^2) f[y]
      - (-2 + 3 b y) (3 b D[f[y], y] + (-2 + 3 b y) D[f[y], {y, 2}]);
bc = DirichletCondition[f[y] == 0, True];
{vals, funcs} =
  NDEigensystem[{lhs, bc}, f[y], {y, -ycutoff, ycutoff}, Nless,
   Method -> {"PDEDiscretization" -> {"FiniteElement", {"MeshOptions"
-> {"MaxCellMeasure" -> 0.01, "MeshOrder" -> 2}}}}];
vals