8th Order Runge-Kutta for Integrating System of ODEs (ODE87)

This solver ODE87 is a realization of explicit Runge-Kutta method.  Integrates a system of ordinary differential equations using 8-7 th order Dorman and Prince formulas. See P.J. Prince & J.R. Dorman (1981) High order embedded Runge-Kutta formulae. J.Comp. Appl. Math., Vol. 7. p.67-75.

This is a 8th-order accurate integrator therefore the local error normally expected is O(h^9).  This requires 13 function evaluations per integration step. Some information about method can be found in Hairer, Norsett and Wanner (1993): Solving Ordinary Differential Equations. Nonstiff Problems.  2nd edition. Springer Series in Comput. Math., vol. 8.

Interface to program based on standart MATLAB ode-suite interface but with some restriction.

• [T,Y] = ODE87(ODEFUN,TSPAN,Y0) with TSPAN = [T0 TFINAL] integrates the system of differential equations y’ = f(t,y) from time T0 to TFINAL with initial conditions Y0. Function ODEFUN(T,Y) must return a column vector corresponding to f(t,y). Each row in the solution array Y corresponds to a time returned in the column vector T.
• [T,Y] = ODE87(ODEFUN,TSPAN,Y0,OPTIONS) solves as above with default integration properties replaced by values in OPTIONS, an argument created with the ODESET function. See ODESET for details. Commonly used options are scalar relative error tolerance ‘RelTol’ (1e-6 by default).
• [T,Y] = ODE87(ODEFUN,TSPAN,Y0,OPTIONS,P1,P2…) passes the additional parameters P1,P2,… to the ODE function as ODEFUN(T,Y,P1,P2…), and to all functions specified in OPTIONS. Use OPTIONS = [] as a place holder if no options are set.

Example on using the solver:

Solve the system y’ = Fun(t,y), using the default relative error tolerance 1e-3 and the default absolute tolerance of 1e-6 for each component, and plots the first component of the solution.  The equations we interested to solve:
y1′ = y(2);
y2′ = -y(1)/sqrt(y(1)^2+y(2)^2)^3;
With initial condtitions y1(0)=.5, y2(0)=0

First form the inline function:
FUN = inline(‘[y(2);-y(1)/sqrt(y(1)^2+y(2)^2)^3]’,’x’,’y’);

Then Call the Solver:
Int = [0 20];   % Interval
IC = [.5 0]’;     % Initial Conditions
[x,y] = ode87(FUN, Int, IC);

Finally visualize the results:
plot(x,y);