d[r_]=r^2 + a^2 -2 M r
rho2[r_,t_]=r^2+a^2 Cos[t]^2
gij=DiagonalMatrix[{rho2[r,t]/d[r], rho2[r,t], Sin[t]^2(r^2+a^2+ 2 M r a^2/rho2[r,t] Sin[t]^2),-(1-2 M r/rho2[r,t])}]
gij[[3,4]]= -2 M r a/rho2[r,t] Sin[t]^2
gij[[4,3]]=gij[[3,4]]

U={r'[s],0,p'[s],ct'[s]}
u=Evaluate[gij /. t->Pi/2].U

out=Solve[{u[[3]]==L,u[[4]]==-e},{p'[s],ct'[s]}]                                                                        

eqp=Collect[First[out][[1]][[2]],{e,L}]
eqt=Collect[First[out][[2]][[2]],{e,L}]

U.Evaluate[gij /. t->Pi/2].U
% /. {ct'[s] -> eqt, p'[s]-> eqp}

Simplify[%]
Solve[%==-1,r'[s]]
r'[s]^2 /. Last[%] 
Together[%]
Collect[%,r,Simplify] 
eqr=%

NOTE:
eqp:
               2 a e M             L (2 M - r)
Out[67]= ------------------- - -------------------
             2            2        2            2
         r (a  - 2 M r + r )   r (a  - 2 M r + r )

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      2      3
              -2 a L M         e (-2 a  M - a  r - r )
Out[68]= ------------------- - -----------------------
             2            2          2            2
         r (a  - 2 M r + r )     r (a  - 2 M r + r )


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:
                                 2      2        2     2
               2   2 (-(a e) + L)  M   a  (-1 + e ) - L    2 M
Out[69]= -1 + e  + ----------------- + ----------------- + ---
                           3                   2            r
                          r                   r

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

near circular orbit:
Solve[{eqr==0,D[eqr,r]==0},{e,L}]
out=% /. {M->1,r->10,a->.98}   
Out[71]= {{e -> -0.962753, L -> 4.24238}, {e -> 0.962753, L -> -4.24238}, 
 
>    {e -> -0.951972, L -> -3.43292}, {e -> 0.951972, L -> 3.43292}}
                                                                            

eqr /. Join[out[[2]],{M->1,a->.98}]
Plot[-%,{r,5,50}]

T=230
all3={eqp, eqt,eqr} /. {M->1,a->.98, L->-4.24238 ,e->0.965,r->r[s]}
r0=10
dr0=Sqrt[eqr] /. {M->1,a->.98,  L->-4.24238 ,e->0.965,r->r0}
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,329.16},AspectRatio->Automatic]

eqr /. Join[out[[2]],{M->1,a->.98}]
Plot[-%,{r,5,50}]

T=450
all3={eqp, eqt,eqr} /. {M->1,a->.98, L->3.43292 ,e->0.955,r->r[s]}
r0=10
dr0=Sqrt[eqr] /. {M->1,a->.98,  L->3.43292,e->0.955,r->r0}
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]
ParametricPlot[Evaluate[{r[s] Cos[p[s]], r[s] Sin[p[s]]} /. Last[solution]],{s,0,10 T},AspectRatio->Automatic]                                         
eqr /. Join[out[[4]],{M->1,a->.98}]
Plot[-%,{r,5,50}]

eqp/eqt /. Join[out[[2]],{M->1,a->.98,r->10}]
Out[100]= -0.0326341
eqp/eqt /. Join[out[[4]],{M->1,a->.98,r->10}]
Out[101]= 0.0306722
