There are lists and descriptions of methods, including implicit one-step methods:

• Under the "Details" of NDSolve

• In the ODE Integration Methods chapter of the NDSolve Overview

LSODA is the default solver for ODE's. IDA for DAEs. Not sure about other forms of DEs.

According to the Overview, "FixedStep" is restricted to one-step methods, so you can't use BDF or LSODA with it. A list of such submethods is in the docs. for NDSolve. Strangely, the FixedStep section gives an example that claims to show how to use Method -> "BDF" using fixed steps, but the steps are in fact of different sizes. You can use BDF for time integration in a PDE, but I do not know how to fix the step size or to restrict the minimum step size.

I do not know of any implementation of Rosebrock methods in Mathematica/WL. How to implement your own methods is described in the Plugin Framework.

Suppose that I don't want to use automatic features of NDSolve, but just prefer fixed, uniform spatial and temporal grids, combined with my choice of the spatial and temporal discretization orders, plus my choice of WorkingPrecision. How to do this? I suspect I might choose SpatialDiscretization->MaxSteps == SpatialDiscretization->MinSteps == (my value) to eliminate automatic spatial grid selection. Is this correct? But how to set the other options so that PrecisionGoal and AccuracyGoal are effectively not used by NDSolve? Given the fragility of NDSolve adaptive calculations, the above approach might be the best one, if only no criptic error messages occur; I could test the accuracy of the procedure by comparing the actual errors with some known solutions, instead of relying on the adaptive algorithm. Leslaw

I don't understand the connection of your remarks about errors and error message to the problem at hand. I get no errors with this:

Clear[foo, sol, x, w, psi];
ClearSystemCache["Numeric"]; (* Why? IDK. Seems to fix a bug when re-run. *)
nstep = 0; neval = 0; (* steps, PDE evaluations *)
PrintTemporary@Dynamic@{Clock@Infinity, foo, nstep, neval}; (* monitors progress *)
xmax = 10;
wmax = 10;
xmin = 10^-20;
psiinit[w_, x_] := (Exp[-w^2]/Sqrt[\[Pi]] - w*Erfc[w])*x;
sol = NDSolve[{2*x*D[psi[w, x], x] == 
    D[psi[w, x], {w, 2}] + 2*w*D[psi[w, x], w] - (x*psi[w, x])^2, 
   psi[w, xmin] == psiinit[w, xmin], 
   Derivative[1, 0][psi][0, x] == -x, psi[wmax, x] == 0}, 
  psi, {w, 0, wmax}, {x, xmin, xmax},
  PrecisionGoal -> 20, AccuracyGoal -> 10,
  WorkingPrecision -> 40,
  StepMonitor :> (++nstep; foo = x),
  EvaluationMonitor :> (++neval),
  Method -> {"MethodOfLines", 
    "SpatialDiscretization" -> {"TensorProductGrid", 
      "MaxStepSize" -> 1/5000, "DifferenceOrder" -> 4}}]

Takes a few seconds to set up the first step. Takes several minutes to get to x > xmin. I have only 32GB RAM and NDSolve[] uses around 38.5GB, so much of the slowness is probably due to swapping.

Have you tried "Pseudospectral"?:

Method -> {"MethodOfLines",
  "SpatialDiscretization" -> {
    "TensorProductGrid",
    (*"MinPoints" -> 65,(* or higher *)*)
    "DifferenceOrder" -> "Pseudospectral"}}

A trial run shows that psi approaches zero very rapidly, going below $10^{-20}$. To get 20 digits of precision, the AccuracyGoal needs to be more than the maximum of $20-\log_{10} | \psi |$ throughout the solution domain. I haven't been able to determine the maximum. To achieve, say, AccuracyGoal -> 60, I need WorkingPrecision to be around 400 or more, or some internal numerical glitch happens and I get an error like the one you described (but did not provide code for!). It also takes a long time. Note: Compilers do not give run-time errors, which is what this is. I used to get an uncaught exception and a "core dump" that could be used with a debugger, if the program had been compiled with debugger metadata. Mathematica gives you a stack trace, but only of user-level expressions. I don't get to trace the numerical glitch, which probably occurs in an untraceable library anyway. It would be nice if it indicated the "non-numerical value" it encountered, instead of just giving the time. It would be a good hint to what went wrong.

Have you tried WorkingPrecision -> 40 or some such setting for it?

Sorry, the computer ate my post. I don't have time to rewrite it. :/

Here is the full precision (53-bits, 15.95 digits):

Style[psi[0, 7] /. sol, PrintPrecision -> 20]
(*  {1.1442768817404132}  *)

Read http://reference.wolfram.com/language/tutorial/Numbers.html for more about kinds of numbers.