Je sais que je me répète. Mais mon message n'était manifestement pas passé.
Qu'un schéma ait un ordre supérieur ne le rend pas forcément meilleur. À autres qualités (stabilité, conservation de l'énergie) équivalentes par rapport à un autre schéma d'ordre inférieur, s'il débouche sur une simulation aussi rapide à pas de temps identique, alors oui, il est clairement mieux. Si le schéma d'ordre inférieur est capable de fournir une simulation de même durée avec un pas de temps nettement inférieur, il se peut très bien que celui-ci soit meilleur.
Pour l'histoire de faire les choses de manière explicite ou implicite, le choix se fait en fonction du problème à modéliser. Par exemple, pour modéliser des phénomènes dominés par l'advection, on préfère des schémas explicites (pour des raisons de précision), alors que pour des phénomènes essentiellement diffusifs, on préfère plutôt des schémas implicites (pour des raisons de stabilité). Pour le cas qui nous intéresse, on a tout intérêt à utiliser un schéma symplectique, i.e. qui trouve le "juste milieu" entre un schéma explicite "pur" et un schéma implicite "pur", sinon une planète va rapidement, et incorrectement, soit se faire éjecter, soit se crasher sur le Soleil.
On a ici une équation très simple. Donc, ça ne me semble pas rentable d'utiliser un intégrateur d'ordre trop élevé, qui rend les calculs lourds et donc empêche de prendre un pas de temps petit. Il ne faut pas non plus avoir peur de la condition CFL : les planètes sont supposées se mouvoir à des vitesses ne variant que très peu, sans pic de vitesse imprévisible, en tout cas avant l'horizon de Lyapunov ; donc, il suffit de choisir un pas de temps suffisamment petit, et tout devrait bien se passer.