The dsolve command with the options numeric and method=classical or method=classical[choice] finds a numerical solution by using one of the classical numerical methods as described in the following.

These methods use a fixed step size, provide no error estimation or correction, and are provided primarily for educational purposes. If practical numerical solution of ODE initial value problems (IVP) is required, it is suggested that any of the other methods mentioned in dsolve[numeric,IVP] help page be used instead.

A number of classical methods are available through use of the choice option in the method=classical[choice] selection. The default is $\mathrm{foreuler}$, the forward Euler method.

The available choices of the classical method are described for the problem $y\'=f\left(t,y\right)$, where $Y_i$ is the estimated value of the solution at the time $t_i$, $h$ is the fixed step size $t_i-t_{i-1}$, and for each formula the value of the solution $Y_{n+1}$ at time $t_{n+1}$ is being computed.

$\mathrm{foreuler}$ is the forward Euler method specified by the equation:

$\mathrm{heunform}$ is the Heun formula (also known as the improved Euler method). It uses the forward Euler method to predict the solution, and then applies the trapezoid rule as a corrector. It is specified by the pair of equations:

$\mathrm{impoly}$ is the improved polygon method (also known as the modified Euler method), as specified by the equation:

$\mathrm{rk2}$ is the second-order classical Runge-Kutta method, as specified by:

Note: This is the same as the Heun formula.

$\mathrm{rk3}$ is the third-order classical Runge-Kutta method, as specified by:

$\mathrm{rk4}$ is the fourth-order classical Runge-Kutta method, as specified by:

This is not to be confused with $\mathrm{method}=\mathrm{rkf45}$, which uses a Fehlberg fourth-fifth order Runge-Kutta method.

$\mathrm{adambash}$ is the Adams-Bashforth method (a predictor method), as specified by:

$\mathrm{abmoulton}$ is the Adams-Bashforth-Moulton method (a predictor-corrector method), as specified by:

where $f\left(t_{n+1}\,Y_{n+1}\right)$ is found by first applying the Adams-Bashforth method (the predictor), then using the above Adams-Bashforth-Moulton method (the corrector).

Note that both $\mathrm{adambash}$ and $\mathrm{abmoulton}$ are multipoint methods that require the initial condition and three other equally spaced starting values. These starting values are obtained by computing the first 3 steps using the $\mathrm{rk4}$ method. In addition, since fixed spacing is required for use of these methods, the final step used in obtaining requested solution values also uses the $\mathrm{rk4}$ method (as there is no guarantee that solution values will only be requested at exact multiples of the step size).

The following options are available for the classical method:

'output'

Specifies the desired output from dsolve, and the known option specifies user-defined known functions. For more information, see dsolve[numeric].

'maxfun'

Specifies a maximum on the number of evaluations of the right-hand side of the first order ODE system. This option can be disabled by specifying $\mathrm{maxfun}=0$. The default value for classical methods is $50000$. For more information, see dsolve[maxfun].

'number', 'procedure', 'start', 'initial', and 'procvars'

These options are used to specify the IVP using procedures. For more information, see dsolve[numeric,IVP].

'startinit','implicit', and 'optimize'

These options control the method and behavior of the computation. For more information on the first two, see dsolve[numeric,IVP], for the last, see dsolve[numeric].

'stepsize'

Specifies the static size of each step. Note that classical does not use error correction estimates to adapt the step size. The default step size is $h=\min \left(\frac{\mathrm{Tf}}{3}-\frac{\mathrm{T0}}{3}\,0.005\right)$ for forward integration and $h=\max \left(\frac{\mathrm{Tf}}{3}-\frac{\mathrm{T0}}{3}\,\mathrm{-0.005}\right)$ for backward integration, respectively, where the problem is solved over the interval from T0 to Tf.

'corrections'

A positive integer that specifies the number of corrections that should be applied for the $\mathrm{abmoulton}$ method at each step. It is recommended that this number not exceed 4 (values greater than 4 generally do not produce more accurate results, and result in longer running times).

Results can be plotted using the odeplot function in the plots package.

$\mathrm{dsys1}≔\mathrm{diff}\left(y\left(t\right)\,\mathrm{`\$`}\left(t\,3\right)\right)=y\left(t\right)+\mathrm{diff}\left(x\left(t\right)\,t\right),\mathrm{diff}\left(x\left(t\right)\,\mathrm{`\$`}\left(t\,2\right)\right)=x\left(t\right)y\left(t\right)-\mathrm{diff}\left(y\left(t\right)\,\mathrm{`\$`}\left(t\,2\right)\right)$

$\mathrm{dsys1}≔\frac{{\ⅆ}^3}{\ⅆt^3}y\left(t\right)=y\left(t\right)+\frac{\ⅆ}{\ⅆt}x\left(t\right),\frac{{\ⅆ}^2}{\ⅆt^2}x\left(t\right)=x\left(t\right)y\left(t\right)-\frac{{\ⅆ}^2}{\ⅆt^2}y\left(t\right)$

$\mathrm{init1}≔y\left(0\right)=1,\mathrm{D}\left(y\right)\left(0\right)=2,{\mathrm{D}}^{\left(2\right)}\left(y\right)\left(0\right)=-1,x\left(0\right)=4,\mathrm{D}\left(x\right)\left(0\right)=4.3$

$\mathrm{init1}≔y\left(0\right)=1,\mathrm{D}\left(y\right)\left(0\right)=2,{\mathrm{D}}^{\left(2\right)}\left(y\right)\left(0\right)=\mathrm{-1},x\left(0\right)=4,\mathrm{D}\left(x\right)\left(0\right)=4.3$

$\mathrm{ans1}≔\mathrm{dsolve}\left(\left\{\mathrm{dsys1}\,\mathrm{init1}\right\}\,\mathrm{numeric}\,\mathrm{method}=\mathrm{classical}\left[\mathrm{abmoulton}\right]\,\mathrm{corrections}=2\right)$

$\mathrm{ans1}≔\mathbf{proc}\left(\mathrm{x\_classical}\right)...\mathbf{end\; proc}$

$\mathrm{ans1}\left(1.0\right)$

$\left[t=1.0\,x\left(t\right)=13.1246093792716\,\frac{\ⅆ}{\ⅆt}x\left(t\right)=18.7876915350014\,y\left(t\right)=3.76283413230601\,\frac{\ⅆ}{\ⅆt}y\left(t\right)=5.31760397608216\,\frac{{\ⅆ}^2}{\ⅆt^2}y\left(t\right)=10.2504683623596\right]$

$\mathrm{deq2}≔\mathrm{diff}\left(y\left(x\right)\,\mathrm{`\$`}\left(x\,2\right)\right)=y\left(x\right)$

$\mathrm{deq2}≔\frac{{\ⅆ}^2}{\ⅆx^2}y\left(x\right)=y\left(x\right)$

$\mathrm{init2}≔y\left(0\right)=1,\mathrm{D}\left(y\right)\left(0\right)=1$

$\mathrm{init2}≔y\left(0\right)=1,\mathrm{D}\left(y\right)\left(0\right)=1$

$\mathrm{Digits}≔20\:$

$\mathrm{ans2}≔\mathrm{dsolve}\left(\left\{\mathrm{deq2}\,\mathrm{init2}\right\}\,\mathrm{numeric}\,\mathrm{method}=\mathrm{classical}\left[\mathrm{heunform}\right]\,\mathrm{output}=\mathrm{Array}\left(\left[0.01\,0.1\,1\right]\right)\,\mathrm{stepsize}=0.001\right)$

$\mathrm{ans2}≔\left[\begin{array}{c}\left[\begin{array}{ccc}x& y\left(x\right)& \frac{\ⅆ}{\ⅆx}y\left(x\right)\end{array}\right]\\ \left[\begin{array}{ccc}0.01& 1.0100501654020131715& 1.0100501654020131715\\ 0.1& 1.1051708996699415880& 1.1051708996699415880\\ 1.& 2.7182813757517608399& 2.7182813757517608399\end{array}\right]\end{array}\right]$

$\mathrm{Digits}≔10\:$

$\mathrm{deq3}≔\left\{\mathrm{diff}\left(y\left(t\right)\,\mathrm{`\$`}\left(t\,4\right)\right)-\mathrm{diff}\left(y\left(t\right)\,\mathrm{`\$`}\left(t\,3\right)\right)=y\left(t\right)t^2\right\}$

$\mathrm{deq3}≔\left\{\frac{{\ⅆ}^4}{\ⅆt^4}y\left(t\right)-\frac{{\ⅆ}^3}{\ⅆt^3}y\left(t\right)=y\left(t\right)t^2\right\}$

$\mathrm{init3}≔\left\{y\left(0\right)=3.56\,\mathrm{D}\left(y\right)\left(0\right)=12\,{\mathrm{D}}^{\left(2\right)}\left(y\right)\left(0\right)=-4\,{\mathrm{D}}^{\left(3\right)}\left(y\right)\left(0\right)=6.544\right\}$

$\mathrm{init3}≔\left\{y\left(0\right)=3.56\,\mathrm{D}\left(y\right)\left(0\right)=12\,{\mathrm{D}}^{\left(2\right)}\left(y\right)\left(0\right)=\mathrm{-4}\,{\mathrm{D}}^{\left(3\right)}\left(y\right)\left(0\right)=6.544\right\}$

$\mathrm{ans3}≔\mathrm{dsolve}\left(\mathrm{deq3}\mathbf{union}\mathrm{init3}\,\mathrm{numeric}\,\mathrm{method}=\mathrm{classical}\left[\mathrm{rk3}\right]\,\mathrm{output}=\mathrm{listprocedure}\right)\:$

$\mathrm{sol3}≔\mathrm{eval}\left(y\left(t\right)\,\mathrm{ans3}\right)$

$\mathrm{sol3}≔\mathbf{proc}\left(t\right)...\mathbf{end\; proc}$

$\mathrm{sol3}\left(0\right)$

$3.56$

$\mathrm{sol3}\left(0.56\right)$

$9.87504792480695$

Boyce, W.E., and DiPrima, R.C. Elementary Differential Equations and Boundary Value Problems. 5th ed. New York: Wiley, 1992.

Conte, S.D., and C. de Boor. Elementary Numerical Analysis, An Algorithmic Approach. McGraw-Hill, 1980.

Fox, L., and Mayers, D.F. Numerical Solution of Ordinary Differential Equations for Scientists and Engineers. Chapman & Hall, 1987.

Lambert, J.D. Computational Methods in Ordinary Differential Equations. New York: Wiley, 1973.