$L_a=\frac{1}{2}m\dot{r}^2+\frac{1}{2}mr^2\dot{\theta}^2+\frac{GmM}{r}: $ Lagrangian.
$\frac{d}{dt}(\frac{\partial L_a}{\partial \dot{\theta}})-\frac{\partial L_a}{\partial \theta}=0\implies mr^2\dot{\theta}=L$. Angular momentum is conserved.
$\frac{d}{dt}(\frac{\partial L_a}{\partial \dot{r}})-\frac{\partial L_a}{\partial r}=0\implies m\ddot{r} -mr\dot{\theta}^2+\frac{GMm}{r^2}=0$.
The above by the Euler-Lagrange Equations.
$m\ddot{r}=\frac{L\dot{\theta}}{r}-\frac{GMm^2\dot{\theta}}{L}$ swapping in $L$.
$m\frac{d}{d\theta}(\dot{r})=\frac{L}{r}-\frac{GMm^2}{L}$ dividing by $\dot{\theta}$
$m\dot{r}=m\dot{\theta}\frac{dr}{d\theta}=\frac{L}{r^2}\frac{dr}{d\theta}=\frac{d}{d\theta}(\frac{-L}{r})$ by the chain rule.
Put it all together:
$\frac{d^2}{d\theta^2}(\frac{1}{r})+\frac{1}{r}=\frac{GMm^2}{L^2}$
So $\frac{1}{r}=c_1\cos{\theta}+c_2\sin{\theta}+\frac{GMm^2}{L^2}$
Let's define $\theta =0$ as the aphelion, $p$.This is the sum of the semi-major axis and the distance from the focus to the center of the ellipse, so $p=a+c$. The perihelion, occurring when $\theta=\pi$, has length $a-c$. Since the radial velocity is $0$ at the aphelion, $c_2=0$. So:
$\frac{1}{a+c}=c_1+\frac{GMm^2}{L^2}$ ,$\frac{1}{a-c}=-c_1+\frac{GMm^2}{L^2}$
Then $\frac{2a}{a^2-c^2}=2\frac{GMm^2}{L^2}$ ,And $\frac{-2c}{a^2-c^2}=2c_1.$
By definition, eccentricity $\epsilon=c/a$.
$\frac{GMm^2}{L^2}=\frac{a}{a^2(1-\epsilon^2)}$
$\frac{1}{r}=\frac{-c}{a^2-c^2}\cos \theta +\frac{a}{a^2-c^2}=\frac{-a\epsilon}{a^2(1-\epsilon^2)}\cos \theta + \frac{a}{a^2(1-\epsilon^2)}$
So $r= \frac{a(1-\epsilon^2)}{1-\epsilon \cos \theta}$.
For a given power of $r$, averaging by $\theta$ is $\langle r^n\rangle_\theta =\frac{1}{2\pi}\int_0^{2\pi} \frac{a^n(1-\epsilon^2)^n d\theta}{(1-\epsilon \cos \theta)^n}$
To average by time, you need the period. $\tau =\int_0^{\tau}1 dt =\int_0^{2\pi} \frac{1}{\dot{\theta}} d\theta=\int_0^{2\pi} \frac{mr^2d\theta}{L}$
Then the time average is $\langle r^n\rangle_t=\frac{1}{\tau}\int_0^{2\pi}\frac{mr^2\cdot r^nd\theta}{L}$
I'm not sure how to do it without the Residue Theorem of Complex Analysis, but it can be shown $\int_0^{2\pi} \frac{d\theta}{(1 - \epsilon \cos \theta)^n}=\frac{2\pi}{(1-\epsilon^2)^{n/2}}$.
From there you need to find $V= -GMm\langle 1/r\rangle_t = \frac{-2\pi}{\tau} GMm \langle mr/L \rangle_\theta$.
