Integrating the Schwarzschild geodesics equations I am trying to make a graph similar to this one 
from the Wikipedia article about Schwarzschild geodesics:

There is an equation like this: $$ \varphi = \int\frac{dr}{r^2\sqrt{\frac{1}{b^2}-(1-\frac{r_s}{r})\frac{1}{r^2}}}.$$
I did not get how they manage to draw this picture from that equation.
Particular questions:


*

*If it is a polar plot, usually it is opposite, like $r(\varphi)$, but not $\varphi(r)$. How can I draw the polar plot if I have $\varphi(r)$ function?

*What is integration interval here? $ \varphi = \int_{?}^{?}...$
 A: To your first question: Using a parametric plot in $r$ or $\phi$ and plotting $\left\{r\cos(\phi),r\sin(\phi)\right\}$. 
To your second question: that is where things become tricky. Depending on what kind of geodesics comes out the radius shrinks or increases and the angle is also kind of problematic for closed orbits. A formulation with proper time as parameter is better suited for implementation.
This Wolfram Demonstration Geodesics in Schwarzschild Space of Niels Walet has a simple implementation of an integration over proper time.
That being said; with a bit of tweaking one can of course implement the differential equation for $\phi$. J. B. Hartle provides a Mathematica Notebook here as supplement to his book Gravity: An Introduction to Einstein's General Relativity.
A: Just to summarize all good comments and answers.
If I use this simple equation with integration interval from 0 to Rmax
$\Large \varphi(r_s,b,R_{max}) = \Huge\int_{0}^{R_{max}}\frac{dr}{r^2\sqrt{\frac{1}{b^2}-(1-\frac{r_s}{r})\frac{1}{r^2}}}$
I can get only red lines from the picture in my original question, when photon falls in the Black Hole, as M. J. Steil told.

In order to get green lines this equation will not work and solution of differential equation is needed. Good example is here provided by M. J. Steil.
From this exampe I could get this picture where photon makes a loop around black hole and goes back.

