This method takes initial \(\vec{r}, \vec{v}\), calculates classical orbit parameters,
increases mean anomaly and performs inverse transformation to get final \(\vec{r}, \vec{v}\)
The logic is based on formulae (4), (6) and (7) from http://dx.doi.org/10.1007/s10569-013-9476-9