Résolution d'équation différentielle sous python

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. . 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é 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 ?