My previous posts didn't seem to solve your problem. This will be no different. My apologies. I have two remarks, which are the sources of two of the problems you have in the code for NDSolve.

1) Only the first option setting is used (in any function). So the second Method setting is ignored. There are reasons that this behavior is convenient for programs, but none of them are relevant here.

2) The basic data computed for each step by the ODE solvers in NDSolve (including those used in the method of lines) consists of, at a minimum, the current time, value(s) of the dependent variable(s), and the value(s) the temporal derivative(s). This last one means that the value of the derivative is computed for the initial condition, which is why you get the 1/0. error. It is also stored in the output (the interpolating functions returned). In fact, this will happen no matter what method you use, because it happens in the initialization phase before any steps are taken. (Try NDSolve`ProcessEquations[<<NDSolve arguments go here>>], and you will see the error.) Furthermore, you should expect a possible problem when (or if) the integration reaches the other singular point y == 1. I don't have a solution, sorry.

3) Maybe a third point. Not all implicit methods avoid evaluating the ODE at the end points of a time step. Many do avoid them, and I don't think it's an issue here.

I hope someone who can help you finds this discussion.

Appendix

The thingsbelow may be read about in the document (Tech Notes) for NDSolve.

The complete solution data structure components and names:

Table[
 component -> 
  NDSolve`SolutionDataComponent[{t, {x}, {u}, {Derivative[0, 1][
      u]}, {a}, {b}, {c}, {d}},
   component],
 {component, {"Time", "Space", "DependentVariables", 
   "TemporalDerivatives", "DiscreteVariables", "IndexedVariables", 
   "Parameters", "SensitivityParameters"}}]
(*
{"Time" -> t,
 "Space" -> {x},
 "DependentVariables" -> {u},
 "TemporalDerivatives" -> {Derivative[0, 1][u]},
 "DiscreteVariables" -> {a},
 "IndexedVariables" -> {b},
 "Parameters" -> {c},
 "SensitivityParameters" -> {d}}
*)

Not all solution variables are used in each problem:

state = First@
   NDSolve`ProcessEquations[{Derivative[2, 0][u][x, t] == 
      Derivative[0, 1][u][x, t], u[x, 0] == Exp[-x^2], u[0, t] == 1, 
     u[1, t] == 1/E}, u, {x, 0, 1}, {t, 0, 1}];
state@"VariableData"
(*  {t, {x}, {u}, {Derivative[0, 1][u]}, {}, {}, {}, {}}  *)

Executing the following for state above will show the actual data, at the beginning and then after a time-stepping through to t == 0.1:

state@"SolutionData"["Forward"]

NDSolve`Iterate[state, 0.1]

state@"SolutionData"["Forward"]

1) Here's an example of the usage of Method adapted from the advice in the documentation for NDSolve. Note that it — um — "runs" but does not really solve your problem. And it only runs until it hits the other singular point y == 1, at which it blows up.

y0 = 0 + 10^-100; (* you wish to avoid this, I know *)
wmax = 10;
ystep = 1/10;
difford = 2;
wprec = 68;
pdesol = 
 NDSolve[{D[psi[w, y], {w, 2}] + 2*w*D[psi[w, y], {w, 1}] - 
     4*y*(1 - y)*D[psi[w, y], {y, 1}] - 4*(psi[w, y]^2 - psi[w, y]) == 0
   , psi[w, y0] == Erf[w]
   , psi[0, y] == 0
   , psi[wmax, y] == Erf[wmax]}
  , psi, {w, 0, wmax}, {y, y0, 1}
  , WorkingPrecision -> wprec, PrecisionGoal -> 10
  , Method -> {"MethodOfLines", 
    "DifferentiateBoundaryConditions" -> True, 
    "SpatialDiscretization" -> {"TensorProductGrid"
      , "Coordinates" -> {N[Subdivide[0, wmax, 20]]}
      , "DifferenceOrder" -> 2}
    , Method -> {"FixedStep", Method -> {"ImplicitRungeKutta"
        , "DifferenceOrder" -> difford
        , "ImplicitSolver" -> {"Newton", AccuracyGoal -> 20, 
          PrecisionGoal -> 20, "IterationSafetyFactor" -> 1}}}}
  , StartingStepSize -> ystep]

2) I know a lot about how NDSolve works, but not everything. Either I don't know how to make it do what you describe, or NDSolve is not set up to do it.

3) Any IRK method that has a leading 0 in the c-vector (of standard Butcher table) evaluates the ode at the beginning of the time step. And any RK method whose a-matrix has nonzero entries on the diagonal or above is called "implicit." It means the algebraic equations of the RK method cannot be solved explicitly for the updates at each stage in terms of earlier stages. E.g. $y_{k+1} = [f(x_k, y_k) + f(x_k+h,y_{k+1})]\,h/2$ is the implicit trapezoidal rule and evaluates the ode at both endpoints of the time step. Or Lobatto methods such as

{a, b, c} = 
 NDSolve`ImplicitRungeKuttaLobattoIIIACoefficients[
  4, MachinePrecision]

The c-vector is {0., 0.5, 1.}.

The w and y steps used are available with

psi["Coordinates"] /. First@pdesol

As for the method:

1) For the Method option, it's always the first occurrence that is used.

2) For the PDE method, it's the method of lines.

3) I assume you're asking about syntax. For the time integration method (the Method suboption of the main Method option), it's the fixed-step strategy with the submethod specified by the Method sub-suboption of the Method suboption, which is "ImplicitRungeKutta". Note there is only one Method option specified for NDSolve, and it has Method options nested within it. The general syntax is

Method -> { name , methodoptions }

Some NDSolve methods are classified as "strategies" in the documentation. These may apply one or more step-integrators in some way. The step-integrators are specified with Method -> meth. Mathematica supplies default ones if no method is specified.

Here is the single method option, with one item per line, indented to reflect nesting:

Method -> {
  "MethodOfLines"
  , "DifferentiateBoundaryConditions" -> True
  , "SpatialDiscretization" -> {
    "TensorProductGrid"
    , MaxPoints -> 20 (* might be better than "Coordinates" *)
    , "DifferenceOrder" -> 2}
  , Method -> { (* time integration method for MethodOfLines *)
      "FixedStep"
      , Method -> { (* submethod for time steps in FixedStep method *)
          "ImplicitRungeKutta"
          , "DifferenceOrder" -> difford
          , "ImplicitSolver" -> {
              "Newton"
              , AccuracyGoal -> 20
              , PrecisionGoal -> 20
              , "IterationSafetyFactor" -> 1}}}}

If you were asking, how to be sure NDSolve did what it was told to do, that's less clear. I assume it does. I have checked sometimes, usually for class demonstrations, that it does what I'm going to tell the class it does. I might have gotten messages saying in effect, "This method failed. Using method <blank> instead." NIntegrate does that, but I can't remember if NDSolve does or not.

The following are probably less helpful. You'll have to read through "Components and Data Structures" and maybe the "Method Plugin Framework" tutorials to understand it.

This sometimes shows the call to initialize a method. For whatever reason, it doesn't show all the methods, either because they are initialized in another way or because they can't be traced.

Trace[
 NDSolve[ <NDSolve args> ],
 _NDSolve`InitializeMethod,
 TraceInternal -> True]

This returns the result of initialization, an NDSolve`StateData object, about which you will have to read.

state = First@NDSolve`ProcessEquations[ <NDSolve args> ]

This will return some information about the method:

NDSolve`Iterate[state, 1] (* you need to advance the integration at least one step *)
state@"MethodData"[]

One technical question: Is there any way to save in a file and/or print the messages from this forum? I cannot see any option for this.

Leslaw

What do you exactly mean by saying that it "runs"? And how do you find that the solution blows up at y=1? When I try to run this code, I get lots of error messages, and although some Interpolating Function object is created, nothing can be plotted (the plotting box is empty).

One more question: I have an independent numerical procedure that calculates the derivative d(psi[w,y])/ dy at y=0 and returns it as a function psi1[w] (this can be done with a pretty good accuracy, but not analytically). The derivative is finite. Is there any way to combine such calculated derivative with the pde to be solved, so that NDSolve does not encounter the 1/0 problem while testing the pde at y=0? I am trying to replace the pde by the statement

D[psi[w,y],{y,1}]== If[y==0, psi1[w], b] 

where "b" is the symbolic formula for the derivative obtained from the pde. However, this approach also produces errors, namely there is a message that 0 is not a valid variable.

Leslaw

Yes, the disappearance of the temporal derivative at y=0 worries me, too. But I have not invented this equation. It is taken from some literature and certainly possesses a solution which was verified numerically, albeit not with a high accuracy that I want to achieve. Leslaw

Thank you, indeed, I had an error in my PDE (one missing "y"), so that the corrected code is:

ClearAll;
y0=0;
wmax=10;
ystep=1/10;
difford=2;
wprec=8;
pdesol=NDSolve[{D[psi[w,y],{w,2}]+2*w*D[psi[w,y],{w,1}]-4*y*(1-y)*D[psi[w,y],{y,1}]-4*y*(psi[w,y]^2-psi[w,y])==0,psi[w,y0]==Erf[w],psi[0,y]==0,psi[wmax,y]==Erf[wmax]},psi,{w,0,wmax},{y,y0,1},WorkingPrecision->wprec ,Method->{"FixedStep", Method->{"ImplicitRungeKutta","DifferenceOrder"->difford}},
StartingStepSize->ystep,
Method->{"MethodOfLines","DifferentiateBoundaryConditions"->true,"SpatialDiscretization"->{"TensorProductGrid", "MaxPoints"->20, "MinPoints"->20,"DifferenceOrder"->2}}];

However, after correcting this mistake, I still get the same error message: Power: Infinite expression 1/0, etc. plus many following messages. As far as I can guess, they arise as a result of dividing the pde by the factor 4y(1-y) standing at the first derivative of psi with respect to y, at y=0. But I don't understand why such a calculation is performed; it should not occur in implicit methods. I don't know whether my specification of "Method" is correct: is it OK to have two different such specifications? Maybe one of them is ignored? How else I could select methods for spatial and temporal discretizations? Leslaw