V (r) = k*r^(n + 1)

v[u_]:=k/u^(n+1)

k=10

l=3

mu=1

n=1

fakev[u_]:=l^2/2/mu*u^2+v[u]

Plot[{v[1/r],l^2/mu/r^2,fakev[1/r]},{r,0,5},PlotRange->{0,40}]

uequi=NSolve[fakev'[uu]==0,uu,Reals]

requi=1/uu/.uequi

fakev[uu]/.uequi

e=15

Plot[2*mu*e-l^2/r^2-2*mu*v[1/r],{r,0,5},PlotRange->{-20,5}]

utp=NSolve[2*mu*e-l^2/r^2-2*mu*v[1/r]==0,r,Reals]

theta[rmax_]:=l*Integrate[1/Sqrt[2*mu*e-l^2*uu^2-2*mu*v[uu]],{uu,1/rmax,1/0.6438864832989055}]

theta[1]*180/Pi

Plot[180/Pi*theta[rm],{rm,0.6438864832989055,1.041830214874270}]

ParametricPlot[{{rm*Cos[theta[rm]],rm*Sin[theta[rm]]},{rm*Cos[Pi-theta[rm]],rm*Sin[Pi-theta[rm]]},{rm*Cos[Pi+theta[rm]],rm*Sin[Pi+theta[rm]]},{rm*Cos[-theta[rm]],rm*Sin[-theta[rm]]}},{rm,0.6438864832989055,1.041830214874270}]