import numpy as N
from scipy.integrate import odeint

def testodeint( f, y0, ts, true_y, myrtol, myatol):
    """integrate dy/dt=f(y,t) and compare with exact solution in ture_y.
    Returns the maximum global relative error and
    maximum global absolute error. """
    
    y = odeint( f, y0, ts, rtol=myrtol, atol=myatol )

    abserror = max(abs(true_y-y[:,0]))
    relerror = max(abs((true_y-y[:,0])/true_y))

    return y, relerror, abserror

def rhs(y,t):
    return -2*y*t

ts = N.arange(0,2,0.1)
y0 = N.exp(ts[0]**2)
true_y = N.exp(-ts**2)

ys=odeint( rhs, y0, ts)

import pylab
for tol,style in zip([1e-6,1e-7],['x','o']):
    rtol=atol=tol
    ys=odeint( rhs, y0, ts, atol=atol, rtol=rtol)

    pylab.subplot(2,1,1)
    pylab.plot(ts,ys[:,0],style,label='result odeint (rtol=%g, atol=%g)' % (rtol,atol))
    pylab.legend()

    pylab.subplot(2,1,2)
    pylab.plot(ts,ys[:,0]-true_y,style,label='error (rtol=%g, atol=%g)' % (rtol,atol))
    pylab.legend()
    pylab.xlabel('t')
pylab.subplot(2,1,1)
pylab.xlabel('t')
pylab.plot(ts,true_y,label='exact solution')
pylab.legend()
pylab.title('Comparing numerical approximations with the exact solution')
pylab.savefig('odeintcompareexact.png')
pylab.savefig('odeintcompareexact.eps')
pylab.show()