Trajectoire d'une particule en python

Bonsoir
Je souhaite obtenir une fonction python qui trace la trajectoire d'une particule de masse m en fonction de ses conditions initiales ( position et vitesse initiale ) et dont le mouvement est donné par les équations suivantes en coordonnées cartésiennes :
$$
\ddot x = -\frac{x}{(x^2+y^2)^{3/2}},
\qquad
\ddot y = -\frac{y}{(x^2+y^2)^{3/2}}.


$$L'algorithme n'a rien de compliqué on résout ces 2 équations différentielles à l'aide la méthode d'Euler puis on utilise matplotlib pour tracer les courbes souhaitées. J'ai codé une fonction euler qui permet de résoudre une équation différentielle d'ordre 1 sur un intervalle $[a,b]$ :

import numpy as np
def euler(f,a,b,h,y_0):
  X = np.arange (a , b+h , h)
  Y = [y_0]
  Y_p = [ f(a , y_0) ]
  for k in range (1 , len(X)-1):
    Y += [ Y_p[k-1]*h + Y[k-1] ]
    Y_p += [ f(X[k] , Y[k]) ]
  Y += [ Y_p[-1]*h+ Y[-1] ]
  return [X , Y]  

Quelqu'un saurait-il comment modifier cette fonction pour qu'elle résolve une équation différentielle d'ordre 2 ?
Merci de votre aide !