I'm working on some trajectories in general relativity and I've worked myself to these 2 equations \begin{equation} \left(\frac{dr}{d\lambda}\right)^2=\frac{l^2}{b^2}+\frac{l^2}{r^2}\left(1-\frac{2M}{r}\right) \end{equation} \begin{equation} \left(\frac{d\phi}{d\lambda}\right)^2=\frac{l^2}{r^4} \end{equation} where $r$ and $\phi$ are polar coordinates, $M,b$ and $l$ are known constants and $\lambda$ is an affine parameter.
My plan is to solve these numerically to get trajectory in polar coordinates. I'm just not really sure how to go about this, since I've never done it before. I thought about dividing one into another which gives me \begin{equation} \left(\frac{dr}{d\phi}\right)^2=r^4\left(\frac{1}{b^2}+\frac{1}{r^2}\left(1-\frac{2M}{r}\right)\right) \end{equation} and then solve this equation using for example the RK4 method (If there is any better method to do this please let me know).
Is this the right way to go about this?