eqp=(2*a*e*M)/(r*(a^2 - 2*M*r + r^2)) - (L*(2*M - r))/(r*(a^2 - 2*M*r + r^2))
eqt=(-2*a*L*M)/(r*(a^2 - 2*M*r + r^2)) - (e*(-2*a^2*M - a^2*r - r^3))/(r*(a^2 - 2*M*r + r^2))
eqr=-1 + e^2 + (2*(-(a*e) + L)^2*M)/r^3 + (a^2*(-1 + e^2) - L^2)/r^2 + (2*M)/r

pick an e,L,a (these picked to be near circular orbit)
all3={eqp, eqt,eqr} /. {M->1,a->.98, L->3.43292 ,e->0.955,r->r[s]}

pick a consisten intial radial velocity & position
r0=10
dr0=Sqrt[eqr] /. {M->1,a->.98,  L->3.43292,e->0.955,r->r0}

pick an integration range
T=450

solve/plot:
solution =NDSolve[{D[r'[s]^2-all3[[3]],s]==0,r'[0]==dr0,p'[s]==all3[[1]],ct'[s]==all3[[2]],r[0]==r0,p[0]==0,ct[0]==0},{r,p,ct},{s,0,10 T}]

ParametricPlot[Evaluate[{r[s] Cos[p[s]], r[s] Sin[p[s]]} /. Last[solution]],{s,0,T},AspectRatio->Automatic]