from math import*
import numpy as np
import matplotlib.pyplot as plt

N=10000
T=20

def F(t,X):
    Z=[0,0]
    Z[0]=X[1]-1/3*X[0]**3+X[0]
    Z[1]=-X[0]
   
    return np.array(Z)


L= np.linspace(0, T, N + 1)
 
S=[]
S1=[]

X=[np.array([0,0.5])]

for i in range(0,N+1):
    
    
    S.append(X[i][0])
    S1.append(X[i][1])
    Y=[]
    Y.append(X[i])
    Y.append(Y[0]+0.5*T/N*F(L[i],Y[0]))
    Y.append(Y[0]+0.5*T/N*F(L[i]+0.5*T/N,Y[1]))
    Y.append(Y[0]+T/N*F(L[i]+0.5*T/N,Y[2]))
    X.append(Y[0]+T/N*(1/6*F(L[i],Y[0])+ 2/6*F(L[i]+0.5*T/N,Y[1])+ 2/6*F(L[i]+0.5*T/N,Y[2]) +1/6*F(L[i]+T/N,Y[3])))
    
A=[]
A1=[]

W=[np.array([1,3])]

for k in range(0,N+1):
    
    
    A.append(W[k][0])
    A1.append(W[k][1])
    Y=[]
    Y.append(W[k])
    Y.append(Y[0]+0.5*T/N*F(L[k],Y[0]))
    Y.append(Y[0]+0.5*T/N*F(L[k]+0.5*T/N,Y[1]))
    Y.append(Y[0]+T/N*F(L[k]+0.5*T/N,Y[2]))
    W.append(Y[0]+T/N*(1/6*F(L[k],Y[0])+ 2/6*F(L[k]+0.5*T/N,Y[1])+ 2/6*F(L[k]+0.5*T/N,Y[2]) +1/6*F(L[k]+T/N,Y[3])))
    

B=[]
B1=[]

W=[np.array([4,4])]

for k in range(0,N+1):
    
    
    B.append(W[k][0])
    B1.append(W[k][1])
    Y=[]
    Y.append(W[k])
    Y.append(Y[0]+0.5*T/N*F(L[k],Y[0]))
    Y.append(Y[0]+0.5*T/N*F(L[k]+0.5*T/N,Y[1]))
    Y.append(Y[0]+T/N*F(L[k]+0.5*T/N,Y[2]))
    W.append(Y[0]+T/N*(1/6*F(L[k],Y[0])+ 2/6*F(L[k]+0.5*T/N,Y[1])+ 2/6*F(L[k]+0.5*T/N,Y[2]) +1/6*F(L[k]+T/N,Y[3])))

    
    
    
plt.plot(L,S, label="CI 1", color="blue")
plt.plot(L,A, label="CI 2", color="red")
plt.plot(L,B, label="CI 3", color="green")


plt.xlabel("temps")
plt.ylabel("x")
plt.title("Battement de coeur")
plt.legend()
plt.grid()