from scitools.std import *

def cn_system(F0,S0,T,dt):
    """Numerical solution with Crank-Nicolson's method of  
       F'(t)=1-S(t), F(0)=F0, 
       S'(t)=F(t)-1, S(0)=S0, 
       for 0<t<T and with time step dt. 
    """ 

    times=arange(0,T+dt,dt)
    N=len(times)

    FS=matrix(zeros([2,N]))    
    FS[:,0] = matrix([[F0],[S0]])

    A=matrix([[1, dt/2],[-dt/2, 1]])
    AI=A.I 
 
    # Compute numerical approximations
    for n in range(N-1):
        b=matrix([[FS[0,n]+dt-dt/2*FS[1,n]], [FS[1,n]-dt+dt/2*FS[0,n]]])
        FS[:,n+1] = AI*b
     
    # Plotting in Phase plane
    plot(array(FS[0,:]),array(FS[1,:]))   
    xlabel('F')
    ylabel('S')
    axis('equal')