from scitools.std import *

def compare(y0,T,N):
    """Comparing explicit and implict schemes for: 
       y'(t)=y**2(t) for 0 < t < T, 
       y(0)=y0 
    """ 

    dt=float(T)/N
    times=linspace(0,T,N+1)

    ex_appr=zeros(N+1) # explicit approximation
    im_appr=zeros(N+1) # implicit approximation
    
    ex_appr[0]=y0
    im_appr[0]=y0

    # Compute numerical approximations
    for n in range(N): 
        ex_appr[n+1]= ex_appr[n]+dt* ex_appr[n]* ex_appr[n]
        im_appr[n+1]= (1-sqrt(1-4*dt*im_appr[n]))/(2*dt) 

    # Exact solution
    exact=exact_solution(y0,times)

    plot(times,ex_appr)
    hold('on')
    plot(times,im_appr)
    plot(times,exact)
    
    xlabel('t')
    ylabel('y')
    legend('explicit','implicit','exact',)
    hold('off')
   

def exact_solution(y0,t):
    """Exact solution of 
       y'(t)=y**2(t) for 0 < t < T, 
       y(0)=y0 
    """
    return y0/(1-y0*t)