gij=DiagonalMatrix[{(1+2/r), r^2, r^2 Sin[t]^2, -(1-2/r)}]
gij[[1,4]]=2/r
gij[[4,1]]=2/r


gijk=ConstantArray[0,{4,4,4}];
Do[ Do[ Do[
gijk[[i,j,k]]=Switch[k,1,D[gij[[i,j]],r],2,D[gij[[i,j]],t],3,D[gij[[i,j]],p],4,0],{k,4}],{j,4}],{i,4}]

gIJ=Inverse[gij]
Glik=ConstantArray[0,{4,4,4}];
Do[ Do[ Do[
sum=Sum[gIJ[[l,j]] (gijk[[i,j,k]]+gijk[[j,k,i]]-gijk[[k,i,j]]),{j,4}];
Glik[[l,i,k]]=FullSimplify[sum/2],{k,4}],{l,4}],{i,4}]

Glik[[1,All,All]]
Glik[[2,All,All]]
Glik[[3,All,All]]
Glik[[4,All,All]]

Glijk=ConstantArray[0,{4,4,4,4}];
Do[ Do[ Do[ Do [
Glijk[[l,i,j,k]]=Switch[k,1,D[Glik[[l,i,j]],r],2,D[Glik[[l,i,j]],t],3,D[Glik[[l,i,j]],p],4,0],{k,4}],{l,4}],{j,4}],{i,4}]

RDabc=ConstantArray[0,{4,4,4,4}];

Do[ Do[ Do[ Do [
sum1=Sum[Glik[[j,a1,b1]] Glik[[d1,c1,j]]-Glik[[j,a1,c1]] Glik[[d1,b1,j]],{j,4}];
RDabc[[d1,a1,b1,c1]]=FullSimplify[sum1+Glijk[[d1,a1,b1,c1]]-Glijk[[d1,a1,c1,b1]] ],{a1,4}],{b1,4}],{c1,4}],{d1,4}]

Rab=ConstantArray[0,{4,4}];
Do[ Do[
Rab[[a1,b1]]=Sum[RDabc[[d1,a1,b1,d1]],{d1,4}],{a1,4}],{b1,4}]
Rab=Simplify[Rab]

Glik=Glik[[All,All,All]] /. {r->r[s], t->Pi/2}

eqr0=r''[s]+ Glik[[1,1,1]] r'[s]^2 + Glik[[1,2,2]] t'[s]^2 + Glik[[1,3,3]] p'[s]^2 + Glik[[1,4,4]] ct'[s]^2+
2 Glik[[1,1,4]] ct'[s] r'[s]
eqr1=r''[s]+ Glik[[1,1,1]] r'[s]^2 +  Glik[[1,4,4]] ct'[s]^2+2 Glik[[1,1,4]] ct'[s] r'[s]

eqt0=t''[s]+ Glik[[4,1,1]] r'[s]^2 + 2 Glik[[2,1,2]] r'[s] t'[s]+Glik[[2,3,3]] p'[s]^2+ Glik[[4,4,4]] ct'[s]^2
#t'[0]=0 eliminates for t=Pi/2

eqp0=p''[s]+ 2 Glik[[3,1,3]] r'[s] p'[s] + 2 Glik[[3,2,3]] t'[s] p'[s]
eqp= % /. {Cot[t]->0}

#p'[0]=0 eliminates for t=Pi/2

eqct=ct''[s]+ Glik[[4,1,1]] r'[s]^2 +2 Glik[[4,1,4]] r'[s] ct'[s] + Glik[[4,4,4]] ct'[s]^2
#p'[0]=0 & t'[0]=0 assumed


solution =NDSolve[{eqr1==0,r'[0]==0,eqct==0,r[0]==10,ct[0]==0,ct'[0]==1},{r,ct},{s,0,50 }]
ParametricPlot[Evaluate[{r[s], ct[s]} /. Last[solution]],{s,0,39.265},AspectRatio->1]

Plot[Evaluate[r[s] /. Last[solution]],{s,0,39.269},AspectRatio->1]
Plot[Evaluate[ct[s] /. Last[solution]],{s,0,39.269},AspectRatio->1]

 solution =NDSolve[{eqr1==0,r'[0]==.25,eqct==0,r[0]==10,ct[0]==0,ct'[0]==1},{r,ct},{s,0,300 }]
ParametricPlot[Evaluate[{r[s], ct[s]} /. Last[solution]],{s,0,228.49},AspectRatio->1]
Plot[Evaluate[r[s] /. Last[solution]],{s,0,228.49},AspectRatio->1]

 solution =NDSolve[{eqr1==0,r'[0]==.3,eqct==0,r[0]==10,ct[0]==0,ct'[0]==1},{r,ct},{s,0,1000 }]
ParametricPlot[Evaluate[{r[s], ct[s]} /. Last[solution]],{s,0,921.8},AspectRatio->1]

 solution =NDSolve[{eqr1==0,r'[0]==.35,eqct==0,r[0]==10,ct[0]==0,ct'[0]==1},{r,ct},{s,0,1000 }]
ParametricPlot[Evaluate[{r[s], ct[s]} /. Last[solution]],{s,0,1000},AspectRatio->1]

 solution =NDSolve[{eqr1==0,r'[0]==0,eqct==0,r[0]==1.95,ct[0]==0,ct'[0]==.5},{r,ct},{s,0,1000 }]
ParametricPlot[Evaluate[{r[s], ct[s]} /. Last[solution]],{s,0,1000},AspectRatio->1]

#inner light cone solutions
g2={{gij[[1,1]],gij[[1,4]]},{gij[[4,1]],gij[[4,4]]}} /. r->1.95
Solve[{{x, 1} . g2 . {x, 1}} == 0, x]
rp0=x /. Last[%]
rp1=x /. First[%%]

solution =
 NDSolve[{eqr1 == 0, r'[0] == rp0, eqct == 0,
   r[0] == 1.95, ct[0] == 0, ct'[0] == 1}, {r, ct}, {s, 0, 500}]
ParametricPlot[
 Evaluate[{r[s], ct[s]} /. Last[solution]], {s, 0, 154.03},
 AspectRatio -> 1]

solution =
 NDSolve[{eqr1 == 0, r'[0] == rp1, eqct == 0,
   r[0] == 1.95, ct[0] == 0, ct'[0] == 1}, {r, ct}, {s, 0, 5}]
ParametricPlot[
 Evaluate[{r[s], ct[s]} /. Last[solution]], {s, 0, 2},
 AspectRatio -> 1]
 Show[%,%%%,PlotRange->All]


r0=300000/1476.628442581272
ctp0=1/Sqrt[1-2/r0]

 solution =NDSolve[{eqr1==0,r'[0]==0,eqct==0,r[0]==r0,ct[0]==0,ct'[0]==ctp0},{r,ct},{s,0,10000 }]
ParametricPlot[Evaluate[{r[s], ct[s]} /. Last[solution]],{s,0,3216},AspectRatio->1]
ParametricPlot[Evaluate[{r[s], ct[s]-2 Log[Abs[r[s]-2]]} /. Last[solution]],{s,0,3216.45},AspectRatio->1]
ParametricPlot[ Evaluate[{r[s], ct[s] - 2 Log[Abs[r[s] - 2]]} /. Last[solution]], {s,3200, 3216.45}, AspectRatio -> 1]
Plot[ Evaluate[r[s] /. Last[solution]],{s,0,3216.45}]
Plot[ Evaluate[-r'[s] /. Last[solution]],{s,0,3216.45}]


