Si tu cherches vraiment une solution physique, comme l'a dit Gabbro (et d'autres) c'est le problème à $N$ corps, qu'on ne sait pas résoudre dans le cas général.
Si tu cherches une solution purement algorithmique, càd si tu cherches à générer des systèmes aléatoirement (que ce soit pour un jeu, ou pour le plaisir, etc… ) je peux te proposer une solution.
Tu modélises ton système comme un arbre. La racine serait l'étoile (le soleil). Les fils immédiats de la racine seraient les corps en rotation autour du soleil. Les fils des fils de la racine seraient les satellites en rotation autour des planètes, elles mêmes en rotation autour du soleil. Tu peux rajouter des niveaux supérieurs si tu le juges nécessaire.
On considère alors qu'un corps ne subit QUE l'influence de ses ancêtres (dans l'arbre). Si on veut plus de précision, on peut également ajouter les frères, les neveux, les arrières grand-oncle, etc… ça devient vite le bordel et c'est pas franchement nécessaire d'aller aussi loin à mon avis.
Tu peux tirer au hasard la masse et la distance du/de la satellite/planète à sa/son planète/étoile, en évitant les distributions uniforme (peu réaliste), et en prenant une gaussienne par exemple.
Tu considères uniquement l'influence du corps attracteur immédiat (le père dans l'arbre) et tu négliges le reste (grands parents).
Tu prendrais comme paramètres de ta Gaussienne les moyennes effectivement observées par les astronomes (vu le nombre de satellites du système solaire, et des exoplanètes déjà trouvées, il doit y avoir assez de données pour commencer à dégager des statistiques, reste à les trouver).
Tu effectues ton tirage aléatoire de distance/masse du corps orbitant en prenant en compte les caractéristiques du corps autour duquel il orbite (pour éviter d'avoir un Jupiter tournant autour de la Lune).
Tu considères que l'orbite est circulaire. Alors (je ne me suis pas soucié des signes et des directions, juste des normes) :
$$ma = \frac{GMm}{r^2}$$
$$m\frac{v^2}{r}=\frac{GMm}{r^2}$$
$$v = \sqrt{\frac{GM}{r}}$$
Je ne suis plus très sûr de mon cours de méca (si un physicien passe, qu'il me corrige !), mais je crois que c'est ça (pour un mouvement circulaire).
Tu peux bien sûr tirer masse/distance au hasard, mais aussi masse/vitesse, ou encore vitesse/distance. Si tu supposes le mouvement circulaire, tu peux toujours déduire la 3ème inconnue des deux premières grâce à la formule que je t'ai donné ci-dessus.
Ainsi, tu tires la masse et la distance de la Lune au hasard. Tu calcuels sa vitesse. Tu lui fais décrire ce mouvement autour de la Terre (avec cette vitesse). De même pour la Terre, connaissant sa masse et sa distance au soleil (tirée aléatoirement) tu calcules sa vitesse. Tu fais tourner la Terre autour du soleil à cette vitesse. Et tu composes la vitesse de la Terre avec celle de la Lune pour que la Lune suive la Terre dans son mouvement tout en tournant autour d'elle.
Cela permet de comprendre un peu mieux le rôle de l'arbre. Il induit une hirarchie entre les astres. Ce modèle est assez imprécis. Surtout si deux gros corps massifs sont sur des orbites trop proches, voire confondues (!), il donne des résultats aberrants.
Donc tu dois identifier les grandeurs caractéristiques de ton problème (ce qu'on physicien peut faire pour le coup), et en déduire pour deux planètes sur des orbites adjacentes (des nœuds frères), quelle distance minimale doit les séparer pour qu'il ne soit pas absurde de négliger l'influence qu'ils ont l'un sur l'autre par rapport à l'influence du soleil (leur nœud père). Il faudra ainsi adapter ton générateur aléatoire pour qu'il prenne ne compte cette nouvelle contrainte. Il faut aussi s'assurer que la masse de la planète soit petite devant celle du soleil, et que la masse du satellite soit petite devant celle de la planète.
C'est de la grosse cuisine imprécise, mais il n'est pas impossible de sortir quelque chose de réaliste et cohérent en ajustant les paramètres du problème.