# WorkingPrecision and NDSolve Method Plugin Framework

I'm studying some new numerical integration algorithms using NDSolve Method Framework to specify them. When I try to set WorkingPrecision option to other value than MachinePrecision I get the following error message:

NDSolve::bdstep: The Step function for Method -> CRK4 returned {0.03529465924388010,0.07025569022824054}, which is not an acceptable form. >>

Here I apply the first CRK4 definition from the official tutorial for the Van der Pol equation:

Needs["DifferentialEquationsNDSolveProblems"];
system =
NDSolveProblem[
{{Y1'[T] == Y2'[T], Y2'[T] == -Y1[T] + 2 (1 - Y1[T]^2) Y2[T]},
{Y1 == 0, Y2 == 2}, {Y1[T], Y2[T]}, {T, 0, 30}, {}, {}, {}}];
vars = system["DependentVariables"];
time = system["TimeData"];
sol1 =
Flatten[vars /. NDSolve[system, Method -> CRK4, WorkingPrecision -> 16]]


I have little experience working with Wolfram Mathematica and I will be very grateful to get a clue how to fix this issue.

Change the definitions so that the step output is a list:

CRK4[]["Step"[rhs_, h_, t_, x_, xp_]] := Module[{k0, k1, k2, k3},
k0 = h xp;
k1 = h rhs[t + h/2, x + k0/2];
k2 = h rhs[t + h/2, x + k1/2];
k3 = h rhs[t + h, x + k2];
{(k0 + 2  k1 + 2  k2 + k3)/6}]        (* <-- List *)

CRK4[___]["StepInput"] = {"F"["T", "X"], "H", "T", "X", "XP"};
CRK4[___]["StepOutput"] = {"XI"};       (* <-- List *)


Then both the MachinePrecision and WorkingPrecision -> 16 code return comparable solutions:

sol1mp = Flatten[vars /. NDSolve[system, Method -> CRK4, StartingStepSize -> 0.01]];
sol1 = Flatten[vars /. NDSolve[system, Method -> CRK4, StartingStepSize -> 0.01,
WorkingPrecision -> 16]];

Plot[sol1mp - sol1 // Evaluate, Evaluate@system["TimeData"], PlotRange -> All] Now, why does it work? Well, the check takes place in internal code that cannot be traced. So I don't know. I don't know if it's a bug or an error in the example. There is quite a bit of documentation on NDSolve[], but there seems to be a lot that is not explained, too. There may be a reason for wanting the two return types.

CRK4[___]["StepOutput"] = "XI";

it seems the output of the "Step" method needs to be a packed array or you'll get the error. That means machine precision only. If the output form is {"XI"}, then NDSolve seems to tolerate any numeric array of the appropriate dimensions.