Résolution d'équation différentielle sous python pour modéliser une trajectoire

Bonjour, je suis actuellement étudiant en prépa MP et je rencontre des problèmes sur mon TIPE. J'ai pour but de modéliser la trajectoire d'astéroïdes sous python. J'ai déjà réussi pour une modélisation assez simple en utilisant la méthode de runge kutta 4. Pièce jointe 619473 Pièce jointe 619474. Voici mon premier code pour vous donner une idée de mon niveau si cela peut vous aider( Je sais que ce n'est pas très compact mais je ne m'intéresse pas à ca pour l'instant).

J'ai voulu alors modéliser mon problème d'une manière plus compliquée en prenant en compte l'action de certaines autres planètes sur mon astéroïde. En réutilisant la méthode de runge kutta comme avant je fais face à un problème. J'ai pris en compte l'action de la Terre, Jupiter et Saturne ensemble. Je veux d'abord pouvoir représenter comme sur la deuxième photo leurs trajectoires en 3D. Mon problème est que dans ce modèle mes équations sont couplées et je n'arrive pas à écrire un programme sur la méthode de runge kutta qui me permet d'obtenir leurs positions comme sur la première photo. Auriez vous une piste pour m'aider ?

Que vous ayez une ou plusieurs planètes ne changent pas la forme général du problème. Simplement il ne faut pas voir cela comme 3 trajectoires couplées, mais comme un seul modèle contenant les 3 trajectoires.
Là vous avez xT,yT,zT, uxT, uyT, uzT (T pour terre je présume). Si je rajoute la lune, et bien vous aurez non plus 6 variables mais 12 : xT,yT,zT, uxT, uyT, uzT, xL,yL,zL, uxL, uyL, uzL. Et vous pourrez exprimer la dérivée de n'importe laquelle de ces variables en fonctions de ces 12 là.