def ode45(f, tspan, y0, tol=1e-6):
    '''
    tt, yy = ode45(f, tspan, y0, tol)
    ODE324 - ODE IVP sovler by Dr Hale and TW324 class 2017
    INPUTS:  A function handle f = @(t,y) ... such that dy/dt = f(t,y)
             The times span tspan(1) < t <tspane(end)
             The initial condition y(tspan(1)) - y0
             The required accuracy, tol.
    OUTPUTS: The selected time steps tt, and solution values yy.
    '''

    # Initial step size
    h = np.sqrt(tol);

    # Initialise solver
    t    = tspan[0];
    tend = tspan[1];
    yn   = y0;
    flag = 0;

    # Storage
    tt = t;
    yy = yn;

    # Loop through time:
    while (t < tend):

        # RK stages
        s1 = f(t, yn)
        s2 = f(t + (0.25)*h, yn + (0.25)*h * s1)
        s3 = f(t + (3./8.)*h, yn + (3./32.)*h*s1 + (9./32.)*h*s2)
        s4 = f(t + (12./13.)* h, yn + (1932./2197)*h*s1 - \
             (7200./2197.)*h*s2 + (7296./2197.)*h*s3)
        s5 = f(t + h, yn  + (439./216.)*h*s1 - 8.0*h*s2 + \
             (3680./513.)*h*s3 - (845./4104.)*h*s4)
        s6 = f(t + .5 * h, yn - (8./27.)*h*s1 + 2.0*h*s2 - \
             (3544./2565.)*h*s3 + (1859.0/4104.0)*h*s4 - (11./40.)*h*s5)
        p  = 5;
        yn1 = yn + h*((16.0/135.0)*s1 + (6656.0/12825.0)*s3 +\
             (28561.0/56430.0)*s4 - (9.0/50.)*s5 + (2./55.)*s6) ;
        # Error estimate
        err = h * abs((1.0/360.0)*s1 - (128.0/4275.0)*s3 - \
              (2197.0/75240.0)*s4 + (1.0/50.0)*s5 + (2.0/55.0)*s6) ;

        # Step size adjust
        hnew = 0.8*(tol*np.min(np.abs(yn1)/err))**(1.0/(p+1.0)) * h;
        hnew = min(hnew, 2.0*h);           # Don't let h grow too fast.

        if (np.abs(err) > tol):            # Error is too large!
            if (flag == 0):                # Try the new time step
                flag = 1;                  # Set flag to remember failure
            else:                          # Time step failed
                hnew = h/2.0;              # Take half of old one
        else:                              # Error is small enough!
            # Store solution and move to next time
            yn = yn1
            t = t + h;
            yy = np.hstack((yy, yn1));
            tt = np.hstack((tt, t));
            # Successful step. Set flag to happy.
            flag = 0;

        # Choose next step (ensure we finish at the end of the interval!)
        h = min(hnew, tend - t);

    # Transpose for convenience
    tt = tt.T
    yy = yy.T
    return (tt, yy)

def ode23(f, tspan, y0, tol=1e-6):
    '''
    tt, yy = ode23(f, tspan, y0, tol)
    ODE324 - ODE IVP sovler by Dr Hale and TW324 class 2017
    INPUTS:  A function handle f = @(t,y) ... such that dy/dt = f(t,y)
             The times span tspan(1) < t <tspane(end)
             The initial condition y(tspan(1)) - y0
             The required accuracy, tol.
    OUTPUTS: The selected time steps tt, and solution values yy.
    '''

    # Initial step size
    h = np.sqrt(tol);

    # Initialise solver
    t    = tspan[0];
    tend = tspan[1];
    yn   = y0;
    flag = 0;

    # Storage
    tt = t;
    yy = yn;

    # Loop through time:
    while (t < tend):

        s1 = f(t, yn)
        s2 = f(t + h, yn + h * s1)
        s3 = f(t + h/2.0, yn + (h/4.0)*(s1 + s2))
        p  = 3
        yn1 = yn + (h/6.0)*(s1 + 4.0*s3 + s2)
        err = h/3.0 * abs(s1 - 2.0*s3 + s2)           # Error estimate

        # Step size adjust
        hnew = 0.8*(tol*np.min(np.abs(yn1)/err))**(1.0/(p+1.0)) * h;
        hnew = min(hnew, 2.0*h);           # Don't let h grow too fast.

        if (np.abs(err) > tol):            # Error is too large!
            if (flag == 0):                # Try the new time step
                flag = 1;                  # Set flag to remember failure
            else:                          # Time step failed
                hnew = h/2.0;              # Take half of old one
        else:                              # Error is small enough!
            # Store solution and move to next time
            yn = yn1
            t = t + h;
            yy = np.hstack((yy, yn1));
            tt = np.hstack((tt, t));
            # Successful step. Set flag to happy.
            flag = 0;
        # Choose next step (ensure we finish at the end of the interval!)
        h = min(hnew, tend - t);
    # Transpose for convenience
    tt = tt.T
    yy = yy.T
    return (tt, yy)

def plot_functions(tt, yy):
    # Plot first component of solution
    plt.subplot(211)
    plt.plot(tt, yy, 'm-', lw=4.0)
    plt.title('First component of solution')
    # Plot time steps as a function of time
    dtt = tt[2:-1] - tt[1:-2]
    plt.subplot(212)
    plt.plot(tt[2:-1], dtt, 'c-', lw=4.0)
    plt.yscale('log')
    plt.ylabel('Time steps')
    plt.title('No. of time steps = %s'%(np.size(tt)-1))
    plt.show()

if __name__ == "__main__":
    import numpy as np
    import numpy.linalg as npla
    import matplotlib.pyplot as plt

    # making use of ODE23
    tt, yy = ode23(lambda t, y: y**2 - y**3, [0, 2.0/0.001], y0=0.001)
    plot_functions(tt, yy)

    # making use of ODE45
    tt, yy = ode45(lambda t, y: y**2 - y**3, [0, 2.0/0.001], y0=0.001)
    plot_functions(tt, yy)
