In[]:=
Clear[r,t,τ,M]
In[]:=
M=500;​​r0=3M;​​uppertime=1996.733812103724;
In[]:=
s=NDSolve[{t'[τ]==Sqrt[1-2M/r0]/(1-2M/r[τ]),r''[τ]-(M/r[τ]^2)/(1-2M/r[τ])*r'[τ]^2+(M/r[τ]^2)*(1-2M/r0)/(1-2M/r[τ])==0,t[0]==0,r[0]==r0,r'[0]==0},{t,r},{τ,0,uppertime}]
Out[]=
tInterpolatingFunction
Domain: {{0.,2000.}}
Output: scalar
,rInterpolatingFunction
Domain: {{0.,2000.}}
Output: scalar

In[]:=
Plot[Evaluate[{t[τ],r[τ]}/.s],{τ,0,uppertime},PlotStyleAutomatic]
Out[]=
500
1000
1500
2000
1000
2000
3000
4000