Yeah, there aren't too many ways of detemining orbital elements from initial conditions so it is very likely that the code implementing that ends up looking the same. The similarity in variable names and naming_conventions made the two functions look rather similar.

What comes to stability: given just a single pair of (position, velocity) vectors, it's hard to make something that is really stable. If you have several pairs, you can use the method of least squares to try to fit an orbit into observations.

Not using too many arcus trig functions (acos, asin), will give a bit better results near zero.

Additionally, not having to solve the Kepler's equation would be an improvement for near-parabolic trajectories. In other words: it might be preferable to store (true anomaly at epoch, timestamp of epoch) rather than time of periapsis passage or mean anomaly at epoch to avoid this. The issue is that E - e*sin(E) loses numerical precision when e is near 1 (parabolic) and E is near zero. Same applies for the hyperbolic equivalent of Kepler's equation.

Anyway, nice to find people with similar interests! Good luck with your projects.