bonjour à tous

J'essaie d'intégrer le système différentiel à deux variables  

y1'(t) = y2(t)

y2'(t) = (1-y1(t)^2)y2(t) - y1(t)

y1(0) = 2.00861...

y2(0) = 0

à l'aide de la méthode de runge mais les resultats obtenue sont archi faux, quelqu'un aurait il une odée d'ou provient l'erreur ? merci

import numpy as np

def runge (t0,tf,fun,x0,N,n):
    T = np.zeros(N+1)
    T[0]= t0
    X = np.zeros([N+1,n])
    X[0,:] = x0
    Pas = (tf-t0)/N
    for k in range (N):
        k1 = fun(X[k][0],X[k][1])
        k2 = fun(X[k][0]+k1[0]*(Pas/2),X[k][1]+ k1[1]*(Pas/2))
        print(k2/2)
        X[k+1] = X[k] + (Pas/2)*(k2)
        T[k+1] = T[k] + Pas
    return(T,X) 
    
def fun(x,y):
    k = y
    l = (1-x**2)*y - x
    return(np.array([k,l]))
t0 = 0
tf = 6.6632868593231301896996820305
N = 10
n = 2
x0 = np.array([2.00861986087484313650940188,0])
runge(t0,tf,fun,x0,N,n)
print(runge(t0,tf,fun,x0,N,n))