A variable time-step, high accuracy ODE solver: ode45.m

Although the code written above will solve systems of ODEs with error at each timestep of $O(h^3)$ (second order accuracy) it may be preferable to use a higher-order method (they are more efficient) and in particular to use an automatically chosen variable timestep. MATLAB has such a routine built in, called ode45.m, which is based on a fourth-order accurate Runge-Kutta method.

$\phantom{0}$Include a call to ode45.m at the end of your script rk2_many.m with the commands:

yin=y(:,1);

[tout,y45] = ode45('rhs2',t,yin);

This will automatically integrate the system of equations whose right hand sides are held in rhs2.m and provide results at the time-points contained in the vector t. The starting values of $y_1$ and $y_2$ are held in the variable yin .

You should compare the output from ode45.m with that from RK2 by using the hold on switch and using the command plot(tout, y45,'k--') to display the results graphically with black dashed lines.

Note 1: the array y45 has $n$ rows and 2 columns, by contrast with y which has 2 rows and $n$ columns. Thus y45(:,1) contains all values of $y_1$ evaluated at the time-points specified in t while y45(:,2) contains all values of $y_2$. The variable tout is identical to t

Note 2: ode45.m outputs results at the points chosen by the user in the variable t , i.e. $h, 2h, \dots$ in our case. However, calculations are actually performed at internal timepoints that are chosen automatically by the routine to give the required accuracy.

Note 3: The built-in default accuracy of ode45.m ( $10^{-3}$ relative error, $10^{-6}$ absolute error) is usually quite adequate for most purposes. However, the inverted pendulum problem is difficult numerically and for some initial conditions the default accuracy gives poor results. The routine can be made more accurate as follows:

options = odeset('RelTol',1.e-6);

[tout,y45] = ode45('rhs',t,yin,options);