NDSolve Method Plugin Framework

Introduction

The control mechanisms set up for (RefLink[NDSolve, paclet:ref/NDSolve]) enable you to define your own numerical integration algorithms and use them as specifications for the (RefLink[Method, paclet:ref/Method]) option of (RefLink[NDSolve, paclet:ref/NDSolve]).

(RefLink[NDSolve, paclet:ref/NDSolve]) accesses its numerical algorithms and the information it needs from them in an object-oriented manner. At each step of a numerical integration, (RefLink[NDSolve, paclet:ref/NDSolve]) keeps the method in a form so that it can keep private data as needed.

The structure for method data used in (RefLink[NDSolve, paclet:ref/NDSolve]).

(RefLink[NDSolve, paclet:ref/NDSolve]) does not access the data associated with an algorithm directly, so you can keep the information needed in any form that is convenient or efficient to use. The algorithm and information that might be saved in its private data are accessed only through method functions of the algorithm object. The key interaction between a method algorithm object and (RefLink[NDSolve, paclet:ref/NDSolve]) is through a "Step" function.

Required "Step" function used by (RefLink[NDSolve, paclet:ref/NDSolve]).

The following is a list of important properties of the method and its "Step" function that should be defined.

Properties that should be defined for using a method algorithm object obj in (RefLink[NDSolve, paclet:ref/NDSolve]).

The table below gives specifications that can be used for step inputs and/or outputs.

Step input and output specifications.

The solution data is the list of values with components relating to different variable types. The solution data in the active direction can be obtained from an NDSolve`StateData object state using state@"SolutionData"["Active"]. Many of the specifications listed above correspond directly to solution data components.

For input, the specification for the function should be given with the arguments it will be given. These arguments should be solution data components. (RefLink[All, paclet:ref/All]) can be used to indicate that all arguments used by (RefLink[NDSolve, paclet:ref/NDSolve]) should be expected.

For output, if the method only returns some of the values in some cases (e.g. the method object does not change every step) then you can use (RefLink[Null, paclet:ref/Null]) for any output except the increment Δx to indicate that it should not be used.

Classical Runge–Kutta

Here is an example of how to set up and access a simple numerical algorithm.

This defines a method function to take a single step toward integrating an ODE using the classical fourth-order Runge–Kutta method. Since the method is so simple, it is not necessary to save any private data:

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]

This specifies for (RefLink[NDSolve, paclet:ref/NDSolve]) the inputs and outputs for the step function. The ___ template is used to make this work for all instances of CRK4 method objects:

CRK4[___]["StepInput"] = {"F"["T","X"],"H","T","X","XP"};
CRK4[___]["StepOutput"] = "XI";

This is an equivalent specification using the full names:

CRK4[___]["StepInput"] = {"Function"["Time","DependentVariables"],"TimeStep","Time","DependentVariables","TemporalDerivatives"};
CRK4[___]["StepOutput"] = "DependentVariablesIncrement";

This defines a method function so that (RefLink[NDSolve, paclet:ref/NDSolve]) can obtain the proper difference order to use for the method:

CRK4[___]["DifferenceOrder"] := 4

This defines a method function for the step mode so that (RefLink[NDSolve, paclet:ref/NDSolve]) will know how to control time steps. This algorithm method does not have any step control, so you define the step mode to be Fixed:

CRK4[___]["StepMode"] := Fixed

This integrates the simple harmonic oscillator equation with fixed step size:

fixed = NDSolve[{x''[t] + x[t] == 0, x[0] == 1, x'[0] == 0},x,{t,0,2 π}, Method->CRK4]

Generally using a fixed step size is less efficient than allowing the step size to vary with the local difficulty of the integration. Modern explicit Runge–Kutta methods (accessed in (RefLink[NDSolve, paclet:ref/NDSolve]) with (RefLink[Method, paclet:ref/Method])->"ExplicitRungeKutta") have a so-called embedded error estimator that makes it possible to very efficiently determine appropriate step sizes. An alternative is to use built-in step controller methods that use extrapolation. The method "DoubleStep" uses an extrapolation based on integrating a time step with a single step of size h and two steps of size h/2. The method "Extrapolation" does a more sophisticated extrapolation and modifies the degree of extrapolation automatically as the integration is performed, but is generally used with base methods of difference orders 1 and 2.

This integrates the simple harmonic oscillator using the classical fourth-order Runge–Kutta method with steps controlled by using the "DoubleStep" method:

dstep = NDSolve[{x''[t] + x[t] == 0, x[0] == 1, x'[0] == 0},x,{t,0,2 π}, Method->{"DoubleStep", Method->CRK4}]

This makes a plot comparing the error in the computed solutions at the step ends. The error for the "DoubleStep" method is shown in blue:

ploterror[{sol_}, opts___] := Module[{
points = x@"Coordinates"[1] /. sol, 
values = x@"ValuesOnGrid" /. sol}, 
ListPlot[Transpose[{points, values - Cos[points]}], opts]
];
Show[
ploterror[fixed, PlotStyle->RGBColor[1,0,1]],
ploterror[dstep, PlotStyle->RGBColor[0,0,1]],PlotRange->All
]

The fixed step size ended up with smaller overall error mostly because the steps are so much smaller; it required more than three times as many steps. For a problem where the local solution structure changes more significantly, the difference can be even greater.

A facility for stiffness detection is described within "DoubleStep Method for NDSolve".

For more sophisticated methods, it may be necessary or more efficient to set up some data for the method to use. When (RefLink[NDSolve, paclet:ref/NDSolve]) uses a particular numerical algorithm for the first time, it calls an initialization function. You can define rules for the initialization that will set up appropriate data for your method.

Initializing a method from (RefLink[NDSolve, paclet:ref/NDSolve]).

NDSolve`InitializeMethod is protected, so to attach rules to it, you would need to unprotect it first. It is better to keep the rules associated with your method. A tidy way to do this is to make the initialization definition using (RefLink[TagSet, paclet:ref/TagSet]) as shown earlier.

As an example, suppose you want to redefine the Runge–Kutta method shown earlier so that instead of using the exact coefficients 2, 1/2, and 1/6, numerical values with the appropriate precision are used instead to make the computation slightly faster.

This defines a method function to take a single step toward integrating an ODE using the classical fourth-order Runge–Kutta method using saved numerical values for the required coefficients:

CRK4[{two_, half_, sixth_}]["Step"[rhs_, h_, t_, x_, xp_]] :=Module[{k0,k1,k2,k3},
k0=h xp;
k1=h rhs[t+half h,x+half k0];
k2=h rhs[t+half h,x+half k1];
k3=h rhs[t+h,x+k2];
sixth (k0+two (k1+k2)+k3)]

This defines a rule that initializes the algorithm object with the data to be used later:

CRK4 /: NDSolve`InitializeMethod[CRK4, stepmode_, sd_, rhs_, state_, opts___] := Module[{prec},
prec = state@"WorkingPrecision";
CRK4[N[{2, 1/2, 1/6}, prec]]]

Saving the numerical values of the numbers gives between 5 and 10 percent speedup for a longer integration using "DoubleStep".

Crank–Nicolson

This is an example of how to set up an implicit method. To solve the system of ODEs x'=f(t,x), the scheme for a time step of size h is (Subscript[x, n+1])= (Subscript[x, n])+(h/2)(f(t, (Subscript[x, n])) + f(t+h, (Subscript[x, n+1]))), where (Subscript[x, n])= x(t) =f(t,x) and (Subscript[x, n])=x(t+h). In general, for nonlinear f, the equations need to be solved with Newton iteration.

Instead of using a true Newton iteration where the Jacobian is updated at each step, the Jacobian will be computed once and its LU decomposition will be used repeatedly to solve the linear equations. As long as time steps are not too large, this will generally converge as well as the full Newton method and each iteration beyond the first is substantially cheaper.

To avoid possible problems with the Newton iteration not converging, the method will be set up to take method options that control the maximum number of iterations and the convergence tolerance.

Define options and method initialization for the Crank–Nicolson method:

Options[CrankNicolson] = {MaxIterations->5, Tolerance->Automatic};

CrankNicolson /: NDSolve`InitializeMethod[CrankNicolson, stepmode_, sd_, rhs_, state_, OptionsPattern[CrankNicolson]] := Module[{prec, rtol, maxit},
maxit = OptionValue[MaxIterations];
prec = state@"WorkingPrecision";
rtol = OptionValue[Tolerance];
If[rtol === Automatic, rtol = 10^(-prec*3/4)];
CrankNicolson[maxit, rtol]]

Define the step function:

CrankNicolson[maxit_, rtol_]["Step"[f_, h_,t0_, x0_, f0_]] :=Module[{J, LU, t1 = t0 + h,x1, f1,residual, err, done = False, tol = rtol, count = 0},
x1 = x0 + h f0;
f1 = f[t1, x1];
x1 = x0 + (h/2)(f0 + f1);
J = f["JacobianMatrix"[t1,x1]];
LU = IdentityMatrix[Length[x1], SparseArray] - (h/2) J;
LU = LinearSolve[LU];
While[(count ≤ maxit) && !done, 
f1 = f[t1,x1];
residual = x1-x0-(h/2)*(f0 + f1);
err = Norm[residual, Infinity];
If[err < tol,
done = True
(* else *),
x1 = x1 - LU[residual];
count++;
]
];
If[count>maxit, Message[CrankNicolson::cvmit, maxit];x1 = $Failed];
{x1, f1}
];

This specifies for (RefLink[NDSolve, paclet:ref/NDSolve]) the inputs, outputs, difference order, and step mode for the Crank–Nicolson method:

CrankNicolson[___]["StepInput"] = {"F"["T","X"],"H","T","X","XP"};
CrankNicolson[___]["StepOutput"] = {"X", "XP"};
CrankNicolson[___]["DifferenceOrder"] := 2;
CrankNicolson[___]["StepMode"] := "Fixed";

Here is how the method can be called with option settings:

NDSolve[{x''[t] + x[t] == 0, x[0] == 1, x'[0] == 0},x,{t,0,2 π}, Method->{"FixedStep","StepSize"->.1, Method->{CrankNicolson, MaxIterations->1, Tolerance->$MachineEpsilon}}]

For a linear problem, the solution is given by the first iteration. However, for a nonlinear problem, convergence may take longer.

This fails because the nonlinear problem cannot be expected to converge in 1 iteration with such large time steps:

NDSolve[{x''[t] + Sin[x[t]] == 0, x[0] == 1, x'[0] == 0},x,{t,0,2 π}, Method->{"FixedStep","StepSize"->.1, Method->{CrankNicolson, MaxIterations->1, Tolerance->$MachineEpsilon}}]

With more iterations the solution is successful:

NDSolve[{x''[t] + Sin[x[t]] == 0, x[0] == 1, x'[0] == 0},x,{t,0,2 π}, Method->{"FixedStep","StepSize"->.1, Method->{CrankNicolson, MaxIterations->5}}]

Since the Crank–Nicolson method is usually considered regarding solving PDEs, here is an example of the method solving the wave equation. Since this is a linear equation, convergence occurs in 1 iteration so the method is quite fast.

This uses the Crank–Nicolson method to solve the wave equation with periodic boundary conditions:

usol = First[u /. NDSolve[{D[u[t,x],t,t]==D[u[t,x],x,x], u[0,x]==Exp[-x^2], Derivative[1,0][u][0,x]==0, u[t,-10]==u[t,10]},u,{t,0,10}, {x,-10,10},
Method->{"DoubleStep", Method->CrankNicolson}]]

A density plot of the solution:

DensityPlot[usol[-t,x],{x,-10,10},{t,-10,0}]

Invoking Another Method

One of the strengths of the (RefLink[NDSolve, paclet:ref/NDSolve]) framework comes from the capability to nest methods in various ways.

Using submethods.

NDSolve`InitializeSubmethod effectively calls NDSolve`InitializeMethod after separating out the method name from any suboptions that are given. If successful, it returns the method object that you can later use in NDSolve`InvokeMethod.

NDSolve`InvokeMethod effectively constructs a "Step" function for the submethod. Since the arguments for the function are specified by the "StepInput" specification for the submethod, the function f given to NDSolve`InvokeMethod should have all the arguments used by (RefLink[NDSolve, paclet:ref/NDSolve]) and can be specified in the caller "StepInput" specification by "Function"[(RefLink[All, paclet:ref/All])]. Note that for built-in methods, there may not be an explicit "Step" function and NDSolve`InvokeMethod handles these using efficient internal code. The result returned by NDSolve`InvokeMethod is always constructed from the submethod's actual output to match a "StepOutput" of {"TimeStep","SolutionData","MethodObject"}.

Here is an example that simply calls a submethod to provide simple monitoring of the steps and difference order.

Define the method options and the initialization for the method:

Options[MonitorMethod] = {Method->Automatic, "MonitorFunction"->Function[{h,sd,mord}, Print[{"H"->h, "SD"->sd,"DifferenceOrder"->mord}]]};

MonitorMethod /: NDSolve`InitializeMethod[MonitorMethod, stepmode_, sd_, rhs_, state_, OptionsPattern[MonitorMethod]] := Module[{submethod, mf},
mf = OptionValue["MonitorFunction"];
submethod = OptionValue[Method];
If[submethod === Automatic, submethod = "StiffnessSwitching"];
submethod = NDSolve`InitializeSubmethod[MonitorMethod, submethod, stepmode, sd, rhs, state];
MonitorMethod[submethod, mf]];

The step function for the method simply invokes the submethod and applies the monitor function to the result and the difference order:

MonitorMethod[submethod_, mf_]["Step"[f_, h_, sd_, state_]] :=
Module[{res = NDSolve`InvokeMethod[submethod, f, h, sd,state]},
mf[h, sd, submethod["DifferenceOrder"]];
If[SameQ[res[[-1]], submethod], res[[-1]] = Null, res[[-1]] = MonitorMethod[res[[-1]], mf]];
res]

Set the inputs and outputs:

MonitorMethod[___]["StepInput"] = {"Function"[All], "H", "SolutionData", "StateData"};
MonitorMethod[___]["StepOutput"] = {"H", "SD", "MethodObject"};

Derive all other properties from the submethod:

MonitorMethod[submethod_,___][prop_] := submethod[prop];

Use the method:

NDSolve[{x''[t]+ x[t]==0,x[0]==1,x'[0]==1},x,{t,0,1}, Method->MonitorMethod, MaxStepFraction->1]

For a more complicated integration the method can be used with (RefLink[Reap, paclet:ref/Reap]) to collect data.

Collect the logarithm of the step size and difference order as a function of time for integrating the Van der Pol oscillator and plot the solution:

{vsol, {data}} = Reap[NDSolve[(({{(Subscript[y, 1])'[t] == (Subscript[y, 2])[t] , (Subscript[y, 1])[0]==2}, {(Subscript[y, 2])'[t]==1000 (1 - (Subscript[y, 1])[t]^2) (Subscript[y, 2])[t] - (Subscript[y, 1])[t], (Subscript[y, 2])[0]==0}})),
{(Subscript[y, 1]), (Subscript[y, 2])},{t,3000}, Method->{MonitorMethod, "MonitorFunction"->Function[{h, sd, mord}, Sow[{NDSolve`SolutionDataComponent[sd, "Time"], RealExponent[h], mord}]]}]];
Plot[(Subscript[y, 1])[t]/. vsol, {t,0,3000}]

Show the logarithm of the step size and the difference order for steps near the sudden change:

Block[{m, δ = 2},
m = t /. FindRoot[(Subscript[y, 1])[t]/. vsol, {t, 0, 1000}];
Show[ListPlot[{data[[All,{1,2}]], data[[All,{1,3}]]}, PlotRange->{{m- δ, m + δ}, All}], Plot[(Subscript[y, 1])[t]/. vsol, {t,m- δ, m + δ}], ImageSize->Medium]]

Adams Methods

In terms of the (RefLink[NDSolve, paclet:ref/NDSolve]) framework, it is not really any more difficult to write an algorithm that controls steps automatically. However, the requirements for estimating error and determining an appropriate step size usually make this much more difficult from both the mathematical and programming standpoints. The following example is a partial adaptation of the Fortran DEABM code of Shampine and Watts to fit into the (RefLink[NDSolve, paclet:ref/NDSolve]) framework. The algorithm adaptively chooses both step size and order based on criteria described in [SG75].

The first stage is to define the coefficients. The integration method uses variable step-size coefficients. Given a sequence of step sizes (Subscript[h, n-k+1]),(Subscript[h, n-k+2]),…,(Subscript[h, n]), where (Subscript[h, n]) is the current step to take, the coefficients for the method with Adams–Bashforth predictor of order k and Adams–Moulton corrector of order k+1, (Subscript[g, j])(n) such that

where the (Subscript[Φ, j])(n) are the divided differences.

This defines a function that computes the coefficients (Subscript[Φ, j]) and (Subscript[β, j]), along with (Subscript[σ, j]), that are used in error estimation. The formulas are from [HNW93] and use essentially the same notation:

AdamsBMCoefficients[hlist_List]:=Module[{k,h,Δh,brat,β,α,σ,c},
k=Length[hlist];
h=Last[hlist];
Δh=Drop[FoldList[Plus,0,Reverse[hlist]],1];
brat=((Drop[Δh,-1])/(Drop[Δh,1]-h));
β=FoldList[Times,1,brat];
α=(h/(Δh));
σ=FoldList[Times,1,α Range[Length[α]]];
c[0]=Table[(1/q),{q,1,k}];
c[1]=Table[(1/(q (q+1))),{q,1,k}];Do[c[j]=Drop[c[j-1],-1]-((Drop[c[j-1],1] h)/(Δh[[j]])),{j,2,k}];{(First[c[#1]]&)/@Range[0,k],β,σ}]

hlist is the list of step sizes (Subscript[h, n-k]),(Subscript[h, n-k+1,])…,(Subscript[h, n]) from past steps. The constant-coefficient Adams coefficients can be computed once, and much more easily. Since the constant step size Adams–Moulton coefficients are used in error prediction for changing the method order, it makes sense to define them once with rules that save the values.

This defines a function that computes and saves the values of the constant step size Adams–Moulton coefficients:

Moulton[0]=1;
Moulton[m_]:=Moulton[m]=-Sum[Moulton[k]/(1+m-k),{k,0,m-1}]

The next stage is to set up a data structure that will keep the necessary information between steps and define how that data should be initialized. The key information that needs to be saved is the list of past step sizes, hlist, and the divided differences, Φ. Since the method does the error estimation, it needs to get the correct norm to use from the (RefLinkPlain[NDSolve`StateData, paclet:tutorial/NDSolveStateData]) object. Some other data such as precision is saved for optimization and convenience.

This defines a rule for initializing the AdamsBM method from (RefLink[NDSolve, paclet:ref/NDSolve]):

AdamsBM/:
NDSolve`InitializeMethod[AdamsBM,{Automatic,__},sd_, rhs_,ndstate_,opts___]:=Module[{prec,norm,hlist,Φ,mord},
mord=MaxDifferenceOrder/.Flatten[{opts,Options[AdamsBM]}];
If[mord ≠ ∞ && !(IntegerQ[mord] && mord > 0),Return[$Failed]];
prec=ndstate["WorkingPrecision"];
norm=ndstate["Norm"];
hlist={};
Φ={ndstate["SolutionDerivativeVector"["Active"]]};AdamsBM[{{hlist,Φ,N[0,prec] Φ[[1]]},{norm,prec,mord,0,True}}]];

hlist is initialized to {} since at initialization time there have been no steps. Φ is initialized to the derivative of the solution vector at the initial condition since the 0(^th) divided difference is just the function value. Note that Φ is a matrix. The third element, which is initialized to a zero vector, is used for determining the best order for the next step. It is effectively an additional divided difference. The use of the other parts of the data is clarified in the definition of the stepping function.

The initialization is also set up to get the value of an option that can be used to limit the maximum order of the method to use. In the code DEABM, this is limited to 12, because this is a practical limit for machine-precision calculations. However, in the Wolfram System, computations can be done in higher precision where higher-order methods may be of significant advantage, so there is no good reason for an absolute limit of this sort. Thus, you set the default of the option to be ∞.

This sets the default for the MaxDifferenceOrder option of the AdamsBM method:

Options[AdamsBM] = {MaxDifferenceOrder->∞};

Before proceeding to the more complicated "Step" method functions, it makes sense to define the simple "StepMode" and "DifferenceOrder" method functions.

This defines the step mode for the AdamsBM method to always be (RefLink[Automatic, paclet:ref/Automatic]). This means that it cannot be called from step controller methods that request fixed step sizes of possibly varying sizes:

AdamsBM[___]["StepMode"]=Automatic;

This defines the difference order for the AdamsBM method. This varies with the number of past values saved:

AdamsBM[data_]["DifferenceOrder"]:=Length[data[[1,2]]];

Finally, here is the definition of the "Step" function. The actual process of taking a step is only a few lines. The rest of the code handles the automatic step size and order selection following very closely the DEABM code of Shampine and Watts.

This defines the "Step" method function for AdamsBM that returns step data according to the templates described earlier:

AdamsBM[data_]["Step"[rhs_,t_,h_,y_,yp_]]:=Module[{prec,norm,hlist,Φ,Φ1,ns,starting,k,zero,g,β,σ,p,f,Δy,normh,ev,err,PE,knew,hnew,temp},{{hlist,Φ,Φ1},{norm,prec,mord,ns,starting}}=data;
(* Norm scaling will be based on current solution y. *)normh=(Abs[h] temp[#1,y]&)/.{temp->norm};
k=Length[Φ];
zero=N[0,prec];
(* Keep track of number of steps at this stepsize h. *)
If[Length[hlist]>0&&Last[hlist]==h,ns++,ns=1];hlist=Join[hlist,{h}];
{g,β,σ}=AdamsBMCoefficients[hlist];
(* Convert Φ to (SuperStar[Φ]) *)
Φ=Φ Reverse[β];
(* PE: Predict and evaluate *)
p=Reverse[Drop[g,-1]].Φ;
f=rhs[h+t,h p+y];
(* Update divided differences *)
Φ=FoldList[Plus,zero Φ1,Φ];
(* Compute scaled error estimate *)
ev=f-Last[Φ];
err=(g[[-2]]-g[[-1]]) normh[ev];
(* First order check: determines if order should be lowered
even in the case of a rejected step *)
knew=OrderCheck[PE,k,Φ,ev,normh,σ];
If[err>1,
(* Rejected step: reduce h by half, make sure starting mode flag is unset and reset Φ to previous values *)
hnew=(h/2);Δy=$Failed;f=None;starting=False;Φ=data[[1,2]],
(* Sucessful step:
CE: Correct and evaluate *)
Δy=h (p+ev Last[g]);
f=rhs[h+t,y+Δy];temp=f-Last[Φ];
(* Update the divided differences *)
Φ=(temp+#1&)/@Φ;
(* Determine best order and stepsize for the next step *)
Φ1=temp-Φ1;knew=ChooseNextOrder[starting,PE,k,knew,Φ1,normh,σ,mord,ns];hnew=ChooseNextStep[PE,knew,h]];
(* Truncate hlist and Φ to the appropriate length for the chosen order. *)
hlist=Take[hlist,1-knew];
If[Length[Φ]>knew,Φ1=Φ[[Length[Φ]-knew]];Φ=Take[Φ,-knew];];
(* Return step data along with updated method data *)
{hnew,Δy,f,AdamsBM[{{hlist,Φ,Φ1},{norm,prec,mord,ns,starting}}]}];

There are a few deviations from DEABM in the code here. The most significant is that coefficients are recomputed at each step, whereas DEABM computes only those that need updating. This modification was made to keep the code simpler, but does incur a clear performance loss, particularly for small to moderately sized systems. A second significant modification is that much of the code for limiting rejected steps is left to (RefLink[NDSolve, paclet:ref/NDSolve]), so there are no checks in this code to see if the step size is too small or the tolerances are too large. The stiffness detection heuristic has also been left out. The order and step-size determination code has been modularized into separate functions.

This defines a function that constructs error estimates (Subscript[PE, j]) for j==k-2, k-1, and k and determines if the order should be lowered or not:

OrderCheck[PE_,k_,Φ_,ev_,normh_,σ_]:=Module[{knew=k},
(Subscript[PE, k])=Abs[σ[[k+1]] Moulton[k] normh[ev]];If[k>1,
(Subscript[PE, k-1])=Abs[σ[[k]] Moulton[k-1] normh[ev+Φ[[2]]]];
If[k>2,
(Subscript[PE, k-2])=Abs[σ[[k-1]] Moulton[k-2] normh[ev+Φ[[3]]]];If[Max[(Subscript[PE, k-1]),(Subscript[PE, k-2])]<(Subscript[PE, k]),knew=k-1]],
If[(Subscript[PE, k-1])<((Subscript[PE, k])/2),knew=k-1];
];
knew
];

This defines a function that determines the best order to use after a successful step:

SetAttributes[ChooseNextOrder,HoldFirst];
ChooseNextOrder[starting_,PE_,k_,knw_,Φ1_,normh_,σ_,mord_,ns_]:=Module[{knew=knw},
starting=starting&&knew≥k&&k<mord;
If[starting,
knew=k+1;(Subscript[PE, k+1])=0,
If[knew≥k&&ns≥k+1,
(Subscript[PE, k+1])=Abs[Moulton[k+1] normh[Φ1]];
If[k>1,
If[(Subscript[PE, k-1])≤Min[(Subscript[PE, k]),(Subscript[PE, k+1])],
knew=k-1,
If[(Subscript[PE, k+1])<(Subscript[PE, k])&&k<mord,knew=k+1]
],
If[(Subscript[PE, k+1])<((Subscript[PE, k])/2),knew=k+1]
];
];
];
knew
];

This defines a function that determines the best step size to use after a successful step of size h:

ChooseNextStep[PE_,k_,h_]:=
If[(Subscript[PE, k])<(2^(-(k+2))),
2 h,
If[(Subscript[PE, k])<(1/2),h,h Max[(1/2),Min[(9/10),((((1/(2 (Subscript[PE, k])))))^(1/(k+1)))]]]
];

Once these definitions are entered, you can access the method in (RefLink[NDSolve, paclet:ref/NDSolve]) by simply using (RefLink[Method, paclet:ref/Method])->AdamsBM.

This solves the harmonic oscillator equation with the Adams method defined earlier:

asol = NDSolve[{x''[t] + x[t] == 0, x[0] == 1, x'[0] == 0},x,{t,0,2 π}, Method->AdamsBM]

This shows the error of the computed solution. It is apparent that the error is kept within reasonable bounds. Note that after the first few points, the step size has been increased:

ploterror[asol]

Where this method has the potential to outperform some of the built-in methods is with high-precision computations with strict tolerances. This is because the built-in methods are adapted from codes with the restriction to order 12.

LorenzEquations = {
{x'[t] == -3 (x[t] - y[t]), x[0] == 0},{y'[t] == -x[t] z[t] + 53/2 x[t] - y[t],y[0] == 1},{z'[t] == x[t] y[t] - z[t], z[0] == 0}};
vars = {x[t],y[t],z[t]};

A lot of time is required for coefficient computation.

Timing[NDSolve[LorenzEquations, vars, {t,0,20}, Method->AdamsBM]]

This is not using as high an order as might be expected.

In any case, about half the time is spent generating coefficients, so to make it better, you need to figure out the coefficient update.

Timing[NDSolve[LorenzEquations, vars, {t,0,20}, Method->AdamsBM, WorkingPrecision->32]]