Integrating the geodesic equation to obtain Newtonian energy of test particle I am studying General Relativity, and have come across a question that I am finding rather intractable:

In Newtonian Theory, the energy equation for a test particle in orbit around a point  mass is: $$\frac{v^{2}}{2} + \frac{\ell^{2}}{2r^{2}} - \frac{GM}{r} = \mathcal{E}$$ Where $r$ is the radius, $v$ is the radial velocity, $\ell$ is the angular momentum per unit mass, $\mathcal{E}$ is the constant energy per unit mass and $-GM/r$ is the gravitational potential. For the Schwarzchild solution show that the integrated geodesic equation may also be written in the form: $$\frac{v_{s}^{2}}{2} + \frac{\ell_{s}^{2}}{2r^{2}} + \Phi_{s}(r) = \mathcal{E}_{s}$$ Where $v_{s} = \mathrm{d}r/\mathrm{d}\tau$, and $\ell_{s}$ and $\mathcal{E}_{s}$ are constants.

I am not really sure what it means by integrating the geodesic equation. My understading is that the geodesic equation is given by:
$$\frac{\mathrm{d}^{2} x^{\lambda}}{\mathrm{d} \tau^{2}} + \Gamma^{\lambda}_{\mu\nu}\frac{\mathrm{d} x^{\mu}}{\mathrm{d} \tau}\frac{\mathrm{d} x^{\nu}}{\mathrm{d} \tau} = 0$$
Where $\Gamma^{\lambda}_{\mu\nu}$ is the affine connection. However, this leads to a system of equations, none of which, when integrated, yield the form given in the question. 
I'm sure that I am missing something fundamental and simple, but I just cannot find it in any of the references or materials that I have been given.

I thought that perhaps integrating the geodesic equation directly would be the way to go, but I do not end up with the correct solution:
$$\int\left(\frac{\mathrm{d}^{2} x^{\lambda}}{\mathrm{d} \tau^{2}} + \Gamma^{\lambda}_{\mu\nu}\frac{\mathrm{d} x^{\mu}}{\mathrm{d} \tau}\frac{\mathrm{d} x^{\nu}}{\mathrm{d} \tau}\right)\:\mathrm{d}\tau = \mathcal{E}_{s}$$
We note that as integration is distributive over addition, we find:
$$\frac{\mathrm{d}x^{\lambda}}{\mathrm{d}\tau} + \int\left(\Gamma^{\lambda}_{\mu\nu}\frac{\mathrm{d} x^{\mu}}{\mathrm{d} \tau}\frac{\mathrm{d} x^{\nu}}{\mathrm{d} \tau}\right)\:\mathrm{d}\tau = \mathcal{E}_{s}$$
We note that if we let $\lambda = r$, then this becomes proportional to $v_{s}$ and thus I am missing a factor of $v_{s}$.

Edit 2: After the very helpful comments below, I went back and did the question again, this time getting stuck at a different point.
So we find that in the $\lambda = t$ case, the only Christoffel symbols $\Gamma_{\mu\nu}^{t}$ that are non-zero are when $\mu = t$, $\nu = r$ or vice versa. Thus we find:
$$\frac{\mathrm{d}t}{\mathrm{d}\tau} + GM\int\frac{\mathrm{d}r}{\mathrm{d}\tau}\frac{\mathrm{d}t}{\mathrm{d}\tau}\frac{\mathrm{d}\tau}{r(r - 2GM)} = \mathcal{E}_{s}$$
We can rewrite the integral in terms of $r$, to give us:
$$\frac{\mathrm{d}t}{\mathrm{d}\tau}\left(1 - \frac{\ln(r)}{2} + \frac{\ln(r - 2GM)}{2}\right) = \mathcal{E}_{s}$$
Which is nowhere near the form that is asked for?
 A: Let $g_{\mu\nu}$ be the metric given by $$\tag{1}\mathrm{d}s^2=-A(r)\mathrm{d}t^2+B(r)\mathrm{d}r^2+r^2\mathrm{d}\Omega^2.$$
The geodesic equation corresponds to the Euler-Lagrange equation of the action $$S[x]=\int\sqrt{-g_{\mu\nu}\dot x^\mu\dot x^\nu}\mathrm{d}\tau$$
Now either (i) vary this functional directly or (ii) calculate the Christoffel symbols of (1) and write down the geodesic equation. Note that $\delta S/\delta x^\mu=0$ is entirely equivalent to $\ddot x^\mu+\Gamma^\mu{}_{\nu\rho}\dot x^\nu\dot x^\rho=0$. Some calculations give 
$$\tag{2}\frac{\delta S}{\delta t}=0\implies\frac{\mathrm{d}}{\mathrm{d}\tau}[A(r)\dot t]=0$$
$$\tag{3}\frac{\delta S}{\delta \phi}=0\implies \frac{\mathrm{d}}{\mathrm{d}\tau}[r^2\sin^2\theta\dot\phi]=0$$ 
$$\tag{4}\frac{\delta S}{\delta\theta}=0\implies \frac{\mathrm{d}}{\mathrm{d}\tau}[r^2\dot\theta]-r^2\sin\theta\cos\theta\dot\phi^2=0$$ 
We can solve (2) and (3) instantly by integration:
$$\tag{5}\dot t=\frac{\epsilon}{A},\quad\dot\phi=\frac{\ell}{r^2\sin^2\theta}$$
where $\epsilon$ and $\ell$ are constants of integration. Additionally, (4) is easily solved by fixing the equatorial plane at $\theta=\pi/2$. We also have the affine parameter equation $$g_{\mu\nu}\dot x^\mu\dot x^\nu=-1=-A(r)\dot t^2+B(r)\dot r^2+r^2\underbrace{\dot\theta^2}_0-r^2\underbrace{\sin^2\theta}_1\dot\phi^2$$
Plugging everything into this equation, we get 
$$\tag{6}\frac{\epsilon^2}{A}-B\dot r^2-\frac{\ell^2}{r^2}=-1$$
Then, use $$A(r)=1-\frac{2GM}{r}=B(r)^{-1}$$
(the Schwarzschild solution) and the desired equation comes out with a little modification of the symbols to fit those in the OP. 
Generally, in the theory of ODEs, such a procedure is called finding a first integral. This is because we integrated two parts of the geodesic equation, namely (2) and (3). It is common to say that the first order equation (6) was obtained by integrating the geodesic equation.
We would like to note that there is a better method for finding first integrals to the geodesic problem. Let $\xi^\mu$ be a Killing vector of the metric $g_{\mu\nu}$ and $\dot x^\mu$ the tangent vector to a geodesic $x^\mu(\tau)$. Then $$g_{\mu\nu}\dot x^\mu\xi^\nu=\text{const.}$$
i.e. it is independent of $\tau$. This is exactly a first integral and may be used to solve the geodesic problem, but there are no actual integrals involved.
