🐞 What bugs you?
I stumbled across this package when I was curious as to whether there was a relativistic equivalent to the REBOUND C/Python package for solving N-body problems. Then I learnt that this package can only predict orbits (find geodesics) for selected metrics for which there's a closed form solution namely Schwarzchild, Kerr and Kerr-Newman, and could not numerically approximate the metric and use this to find geodesics, hence as all these metrics are for single large bodies deforming spacetime, this package is unsuitable for solving N-body problems.

Is there any intention to one day perhaps change this and make this package capable of numerically approximating metrics for such problems? I realize it would not be easy, you'd have to apply something like the finite element method to numerically approximate the metric and it's likely going to be quite a mammoth undertaking, but is there any interest in attempting this?

🎯 Goal
Essentially I want this package to be more general in its ability to predict the trajectory of celestial bodies, able to be applied even when analytical solutions are not available.

💡 Possible solutions
Use some numerical scheme to solve the Einstein field equations (EFEs), after the use specifies what bodies (e.g. planets, stars, black holes) are on the stage of spacetime. Finite element method is one option.

Then use this metric to solve the geodesic equation to predict trajectories.

Conclusion
I did delete the last part of this feature request template because to tell you the truth, I couldn't implement this myself. I know how to numerically integrate ODEs, but not PDEs, especially not systems of PDEs like the EFEs are. I started this issue mostly just to see if there was any interest in implementing such a feature, and whether it even seemed feasible. After all, this project is written in Python, and Python may be too slow to numerically solve the EFEs in real-time.