Bulirsch-Stoer

% [z, info] = BulirschStoer(dynFun,t,z0,tol)
%
% Solves an initial value problem using the Bulirsch-Stoer method. This
% method is ideal for high-accuracy solutions to smooth initial value
% problems.
%
% Computes z(t) such that dz/dt = dynFun(t,z), starting from the initial
% state z0. The solution at the grid-points will be accurate to within tol.
%
% If the provided grid is insufficient, this function will automatically
% introduce intermediate grid points to achieve the required accuracy.
%
% INPUTS:
% dynFun = function handle for the system dynamics
% dz = dynFun(t,z)
% t = scalar time
% z = [nz,1] = state as column vector
% dz = [nz,1] = derivative of state as column vector
% t = [1,nt] = time grid-point vector
% z0 = [nz,1] = initial state vector
% tol = [nz,1] = error tolerance along each dimension. If tol is a
% scalar, then all dimensions will satisfy that error tolerance.
%
% OUTPUTS: (nt = n+1)
% z = [nz,nt] = solution to the initial value problem
% info
% .error = [nx,nt] = error estimate at each grid point
% .nFunEval = [1,nt] = number of function evaluations for each point
%
% NOTES:
% Implementation details:
% http://web.mit.edu/ehliu/Public/Spring2006/18.304/implementation_bulirsch_stoer.pdf
%