(* kerr.newman.yukterez.net || geodesics.yukterez.net || Simon Tyran, Vienna *)

mt1=Automatic;
mt2={"StiffnessSwitching",Method->{"ExplicitRungeKutta",Automatic}};
mt0=mt2; (* Integrationsmethode *)
wp=32; (* Präzision *)

a=9/10; ℧=4/10;  (* Schwarzlochparameter *)
r0=3; θ0=80π/180; φ0=0; (* Startposition *)
μ=-1; q=0;  (* Partikelparameter *)
v1=0;vrain=v2=0; (* Relativgeschwindigkeit *)
vr0=-(((v1-v2))/(1-v1 v2)); vθ0=0; vφ0=0; v0=(Sqrt[(vr0^2)+((vθ0)^2)+((vφ0)^2)]);  (* Startgeschwindigkeit *)
T=10; (* Integrationszeit *)

rA=N[{1-(Sqrt[1-(a^2)-(℧^2)]), 1+(Sqrt[1-(a^2)-(℧^2)])}]; Horizons->rA
rE=N[{1-(Sqrt[1-(a^2) ((Cos[θ0])^2)-(℧^2)]), 1+(Sqrt[1-(a^2) ((Cos[θ0])^2)-(℧^2)])}]; Ergospheres->rE

initcon=Solve[vr0==(((Sqrt[2]) (Sqrt[1-(v0^2) (μ^2)]) (((q r0 (v0^2) (μ^2) ℧ (Sqrt[2 r0-(℧^2)]))/((Sqrt[(a^2)+(r0^2)]) ((a^2)+(-2+r0) r0+(℧^2))))+((dr0 ((a^2)+2 (r0^2)+(a^2) Cos[2 θ0]))/(2 ((a^2)+(r0^2))))+(((Sqrt[2 r0-(℧^2)]) (dt0-a dφ0 ((Sin[θ0])^2)))/(Sqrt[(a^2)+(r0^2)]))))/(Sqrt[(((a^2)+2 (r0^2)+(a^2) Cos[2 θ0])/((a^2)+(r0^2)))]))&&vθ0==dθ0 (Sqrt[1-(v0^2) (μ^2)]) (Sqrt[(r0^2)+(a^2) ((Cos[θ0])^2)])&&vφ0==(((Sqrt[1-(v0^2) (μ^2)]) ((Sin[θ0])^2) (a q r0 (Sqrt[(a^2)+(r0^2)]) (v0^2) (μ^2) ℧-(1/2) a dr0 (Sqrt[2 r0-(℧^2)]) ((a^2)+2 (r0^2)+(a^2) Cos[2 θ0])+(Sqrt[(a^2)+(r0^2)]) (a dt0 (-2 r0+(℧^2))+dφ0 (((((a^2)+(r0^2)))^2)-(a^2) ((a^2)-2 r0+(r0^2)+(℧^2)) ((Sin[θ0])^2)))))/((Sqrt[(a^2)+(r0^2)]) ((r0^2)+(a^2) ((Cos[θ0])^2)) (Sqrt[((((Sin[θ0])^2) (((((a^2)+(r0^2)))^2)-(a^2) ((a^2)-2 r0+(r0^2)+(℧^2)) ((Sin[θ0])^2)))/((r0^2)+(a^2) ((Cos[θ0])^2)))])))&&-μ==-((2 dr0 dt0 (Sqrt[2 r0-(℧^2)]))/(Sqrt[(a^2)+(r0^2)]))+((dθ0)^2) (-(r0^2)-(a^2) ((Cos[θ0])^2))-(((dr0^2) ((a^2)+2 (r0^2)+(a^2) Cos[2 θ0]))/(2 ((a^2)+(r0^2))))+(dt0^2) (1+((-4 r0+2 (℧^2))/((a^2)+2 (r0^2)+(a^2) Cos[2 θ0])))+((2 a dr0 dφ0 (Sqrt[2 r0-(℧^2)]) ((Sin[θ0])^2))/(Sqrt[(a^2)+(r0^2)]))+((2 a dt0 dφ0 (2 r0-(℧^2)) ((Sin[θ0])^2))/((r0^2)+(a^2) ((Cos[θ0])^2)))+((((dφ0)^2) (-((((a^2)+(r0^2)))^2) ((Sin[θ0])^2)+(a^2) ((a^2)+(-2+r0) r0+(℧^2)) ((Sin[θ0])^4)))/((r0^2)+(a^2) ((Cos[θ0])^2)))&&dt0>0,{dt0,dr0,dθ0,dφ0}] (* Anfangsbedingungssolver *)

DGL={
t''[τ]==(1/(8 ((((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)))^3)))(8 q ℧ (-(a^4) ((Cos[θ[τ]])^4)+((r[τ])^4)) r'[τ]+((8 ((a^2) ((Cos[θ[τ]])^2)+(℧^2) r[τ]-((r[τ])^2)) ((((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)))^2) ((r'[τ])^2))/((Sqrt[-(℧^2)+2 r[τ]]) (Sqrt[(a^2)+((r[τ])^2)])))+8 q ℧ (Sqrt[-(℧^2)+2 r[τ]]) (Sqrt[(a^2)+((r[τ])^2)]) (-(a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)) t'[τ]+16 ((a^2) ((Cos[θ[τ]])^2)+(℧^2) r[τ]-((r[τ])^2)) ((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)) r'[τ] t'[τ]+8 (Sqrt[-(℧^2)+2 r[τ]]) ((a^2) ((Cos[θ[τ]])^2)+(℧^2) r[τ]-((r[τ])^2)) (Sqrt[(a^2)+((r[τ])^2)]) ((t'[τ])^2)-8 (a^2) q ℧ r[τ] ((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)) Sin[2 θ[τ]] θ'[τ]+((8 (a^2) (Sqrt[-(℧^2)+2 r[τ]]) ((((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)))^2) Sin[2 θ[τ]] r'[τ] θ'[τ])/(Sqrt[(a^2)+((r[τ])^2)]))-8 (a^2) ((℧^2)-2 r[τ]) ((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)) Sin[2 θ[τ]] t'[τ] θ'[τ]+8 r[τ] (Sqrt[-(℧^2)+2 r[τ]]) (Sqrt[(a^2)+((r[τ])^2)]) ((((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)))^2) ((θ'[τ])^2)-8 a q ℧ (Sqrt[-(℧^2)+2 r[τ]]) (Sqrt[(a^2)+((r[τ])^2)]) (-(a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)) ((Sin[θ[τ]])^2) φ'[τ]-16 a ((a^2) ((Cos[θ[τ]])^2)+((℧^2)-r[τ]) r[τ]) ((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)) ((Sin[θ[τ]])^2) r'[τ] φ'[τ]-16 a (Sqrt[-(℧^2)+2 r[τ]]) ((a^2) ((Cos[θ[τ]])^2)+(℧^2) r[τ]-((r[τ])^2)) (Sqrt[(a^2)+((r[τ])^2)]) ((Sin[θ[τ]])^2) t'[τ] φ'[τ]+16 (a^3) Cos[θ[τ]] ((℧^2)-2 r[τ]) ((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)) ((Sin[θ[τ]])^3) θ'[τ] φ'[τ]+(Sqrt[-(℧^2)+2 r[τ]]) (Sqrt[(a^2)+((r[τ])^2)]) ((a^4)+(a^4) Cos[4 θ[τ]] (-1+r[τ])+3 (a^4) r[τ]+4 (a^2) (℧^2) r[τ]-4 (a^2) ((r[τ])^2)+8 (a^2) ((r[τ])^3)+8 ((r[τ])^5)+4 (a^2) Cos[2 θ[τ]] r[τ] ((a^2)-(℧^2)+r[τ]+2 ((r[τ])^2))) ((Sin[θ[τ]])^2) ((φ'[τ])^2)),
t'[0]==dt0/.initcon[[1]],
t[0]==0,

r''[τ]==(1/(8 ((((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)))^3)))(-8 q ℧ (Sqrt[-(℧^2)+2 r[τ]]) (Sqrt[(a^2)+((r[τ])^2)]) (-(a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)) r'[τ]+((8 (a^2) q ℧ (Sqrt[-(℧^2)+2 r[τ]]) (-(a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)) ((Sin[θ[τ]])^2) r'[τ])/(Sqrt[(a^2)+((r[τ])^2)]))+((4 ((((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)))^2) ((a^2) Cos[2 θ[τ]] (-1+r[τ])-(a^2) (1+r[τ])+2 r[τ] (-(℧^2)+r[τ])) ((r'[τ])^2))/((a^2)+((r[τ])^2)))-4 q ℧ (2 (a^2) ((Cos[θ[τ]])^2)-2 ((r[τ])^2)) ((a^2)+((r[τ])^2)) (1+(((℧^2)-2 r[τ])/((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)))) t'[τ]+((8 (a^2) q ℧ ((℧^2)-2 r[τ]) ((a^2) ((Cos[θ[τ]])^2)-((r[τ])^2)) ((Sin[θ[τ]])^2) t'[τ])/((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)))-((16 (Sqrt[-(℧^2)+2 r[τ]]) ((a^2) ((Cos[θ[τ]])^2)+(℧^2) r[τ]-((r[τ])^2)) ((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)) r'[τ] t'[τ])/(Sqrt[(a^2)+((r[τ])^2)]))+8 ((a^2) ((Cos[θ[τ]])^2)+(℧^2) r[τ]-((r[τ])^2)) ((a^2)+(℧^2)-2 r[τ]+((r[τ])^2)) ((t'[τ])^2)+8 (a^2) ((((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)))^2) Sin[2 θ[τ]] r'[τ] θ'[τ]+8 r[τ] ((((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)))^2) ((a^2)+(℧^2)-2 r[τ]+((r[τ])^2)) ((θ'[τ])^2)+((8 a q ℧ ((℧^2)-2 r[τ]) ((a^2) ((Cos[θ[τ]])^2)-((r[τ])^2)) ((a^2)+((r[τ])^2)) ((Sin[θ[τ]])^2) φ'[τ])/((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)))-((4 a q ℧ (2 (a^2) ((Cos[θ[τ]])^2)-2 ((r[τ])^2)) (-((((a^2)+((r[τ])^2)))^2) ((Sin[θ[τ]])^2)+(a^2) ((a^2)+(℧^2)+(-2+r[τ]) r[τ]) ((Sin[θ[τ]])^4)) φ'[τ])/((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)))-((8 a (Sqrt[-(℧^2)+2 r[τ]]) ((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)) ((a^2) (-1+r[τ])+(a^2) Cos[2 θ[τ]] (-1+r[τ])+2 r[τ] (-(℧^2)+r[τ]+((r[τ])^2))) ((Sin[θ[τ]])^2) r'[τ] φ'[τ])/(Sqrt[(a^2)+((r[τ])^2)]))-16 a ((a^2) ((Cos[θ[τ]])^2)+((℧^2)-r[τ]) r[τ]) ((a^2)+(℧^2)-2 r[τ]+((r[τ])^2)) ((Sin[θ[τ]])^2) t'[τ] φ'[τ]+((a^2)+(℧^2)-2 r[τ]+((r[τ])^2)) ((a^4)+(a^4) Cos[4 θ[τ]] (-1+r[τ])+3 (a^4) r[τ]+4 (a^2) (℧^2) r[τ]-4 (a^2) ((r[τ])^2)+8 (a^2) ((r[τ])^3)+8 ((r[τ])^5)+4 (a^2) Cos[2 θ[τ]] r[τ] ((a^2)-(℧^2)+r[τ]+2 ((r[τ])^2))) ((Sin[θ[τ]])^2) ((φ'[τ])^2)),
r'[0]==dr0/.initcon[[1]],
r[0]==r0,
θ''[τ]==-(1/(2 ((((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)))^4)))(((2 (a^2) Cos[θ[τ]] ((((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)))^3) Sin[θ[τ]] ((r'[τ])^2))/((a^2)+((r[τ])^2)))-2 (a^2) q ℧ ((℧^2)-2 r[τ]) r[τ] Sin[2 θ[τ]] t'[τ]+(a^2) q ℧ r[τ] ((a^2)+2 (℧^2)+(a^2) Cos[2 θ[τ]]-4 r[τ]+2 ((r[τ])^2)) Sin[2 θ[τ]] t'[τ]+2 (a^2) Cos[θ[τ]] ((℧^2)-2 r[τ]) ((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)) Sin[θ[τ]] ((t'[τ])^2)+4 r[τ] ((((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)))^3) r'[τ] θ'[τ]-2 (a^2) Cos[θ[τ]] ((((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)))^3) Sin[θ[τ]] ((θ'[τ])^2)-4 (a^3) q ℧ Cos[θ[τ]] ((℧^2)-2 r[τ]) r[τ] ((Sin[θ[τ]])^3) φ'[τ]+4 a q ℧ Cot[θ[τ]] r[τ] (-((((a^2)+((r[τ])^2)))^2) ((Sin[θ[τ]])^2)+(a^2) ((a^2)+(℧^2)+(-2+r[τ]) r[τ]) ((Sin[θ[τ]])^4)) φ'[τ]+((2 a (Sqrt[-(℧^2)+2 r[τ]]) ((((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)))^3) Sin[2 θ[τ]] r'[τ] φ'[τ])/(Sqrt[(a^2)+((r[τ])^2)]))-2 a ((℧^2)-2 r[τ]) ((a^2)+((r[τ])^2)) ((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)) Sin[2 θ[τ]] t'[τ] φ'[τ]+((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)) (2 (a^2) Cos[θ[τ]] ((Sin[θ[τ]])^3) (-((((a^2)+((r[τ])^2)))^2)+(a^2) ((a^2)+(℧^2)+(-2+r[τ]) r[τ]) ((Sin[θ[τ]])^2))+((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)) (4 (a^2) Cos[θ[τ]] ((a^2)+(℧^2)+(-2+r[τ]) r[τ]) ((Sin[θ[τ]])^3)-((((a^2)+((r[τ])^2)))^2) Sin[2 θ[τ]])) ((φ'[τ])^2)),
θ'[0]==dθ0/.initcon[[1]],
θ[0]==θ0,
φ''[τ]==(1/(4 ((((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)))^3)))(((4 a q ℧ (-(a^4) ((Cos[θ[τ]])^4)+((r[τ])^4)) r'[τ])/((a^2)+((r[τ])^2)))+((4 a ((a^2) ((Cos[θ[τ]])^2)+(℧^2) r[τ]-((r[τ])^2)) ((((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)))^2) ((r'[τ])^2))/((Sqrt[-(℧^2)+2 r[τ]]) ((((a^2)+((r[τ])^2)))^(3/2))))+((4 a q ℧ (Sqrt[-(℧^2)+2 r[τ]]) (-(a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)) t'[τ])/(Sqrt[(a^2)+((r[τ])^2)]))+((8 a ((a^2) ((Cos[θ[τ]])^2)+((℧^2)-r[τ]) r[τ]) ((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)) r'[τ] t'[τ])/((a^2)+((r[τ])^2)))+((4 a (Sqrt[-(℧^2)+2 r[τ]]) ((a^2) ((Cos[θ[τ]])^2)+(℧^2) r[τ]-((r[τ])^2)) ((t'[τ])^2))/(Sqrt[(a^2)+((r[τ])^2)]))-8 a q ℧ Cot[θ[τ]] r[τ] ((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)) θ'[τ]+((8 a Cot[θ[τ]] (Sqrt[-(℧^2)+2 r[τ]]) ((((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)))^2) r'[τ] θ'[τ])/(Sqrt[(a^2)+((r[τ])^2)]))-8 a Cot[θ[τ]] ((℧^2)-2 r[τ]) ((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)) t'[τ] θ'[τ]+((4 a r[τ] (Sqrt[-(℧^2)+2 r[τ]]) ((((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)))^2) ((θ'[τ])^2))/(Sqrt[(a^2)+((r[τ])^2)]))-((4 (a^2) q ℧ (Sqrt[-(℧^2)+2 r[τ]]) (-(a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)) ((Sin[θ[τ]])^2) φ'[τ])/(Sqrt[(a^2)+((r[τ])^2)]))+((8 ((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)) ((a^4) ((Cos[θ[τ]])^2) (-1+r[τ])-r[τ] ((a^2) ((a^2)+(℧^2)-r[τ])+2 (a^2) ((Cot[θ[τ]])^2) ((a^2)+((r[τ])^2))+((Csc[θ[τ]])^2) (-(a^4)+((r[τ])^4)))) ((Sin[θ[τ]])^2) r'[τ] φ'[τ])/((a^2)+((r[τ])^2)))-((8 (a^2) (Sqrt[-(℧^2)+2 r[τ]]) ((a^2) ((Cos[θ[τ]])^2)+(℧^2) r[τ]-((r[τ])^2)) ((Sin[θ[τ]])^2) t'[τ] φ'[τ])/(Sqrt[(a^2)+((r[τ])^2)]))-Cot[θ[τ]] ((a^2) ((Cos[θ[τ]])^2)+((r[τ])^2)) ((a^2) (3 (a^2)-4 (℧^2)+4 ((a^2)+(℧^2)) Cos[2 θ[τ]]+(a^2) Cos[4 θ[τ]])+16 (a^2) ((Cos[θ[τ]])^2) ((r[τ])^2)+8 ((r[τ])^4)+16 (a^2) r[τ] ((Sin[θ[τ]])^2)) θ'[τ] φ'[τ]+((4 a (Sqrt[-(℧^2)+2 r[τ]]) ((Sin[θ[τ]])^2) (r[τ] (-(a^4)+((r[τ])^4)+(a^2) ((a^2)+(℧^2)-r[τ]) ((Sin[θ[τ]])^2))+((Cos[θ[τ]])^2) (2 (a^2) r[τ] ((a^2)+((r[τ])^2))-(a^4) (-1+r[τ]) ((Sin[θ[τ]])^2))) ((φ'[τ])^2))/(Sqrt[(a^2)+((r[τ])^2)]))),
φ'[0]==dφ0/.initcon[[1]],
φ[0]==φ0
}; (* Differentialgleichung *)

sol=NDSolve[DGL,{t,r,θ,φ},{τ,0,T},
MaxSteps->Infinity,InterpolationOrder->All, WorkingPrecision->wp,Method->mt0,StepMonitor:>(laststep=plunge;plunge=τ; stepsize=plunge-laststep;),Method->{"EventLocator","Event":>(If[stepsize<1*^-4,0,1])}]; (* Integrator *)

ct[τ_]:=Evaluate[t[τ]/.sol][[1]]; (* Koordinatenzeit *)
cr[τ_]:=Evaluate[r[τ]/.sol][[1]]; (* Radialkoordinate *)
cθ[τ_]:=Evaluate[θ[τ]/.sol][[1]]; (* Längengradkoordinate *)
cφ[τ_]:=Evaluate[φ[τ]/.sol][[1]]; (* Breitengradkoordinate *)

vf[τ_]:=NSolve[vR==(((Sqrt[2]) (Sqrt[1-(μ^2) (vT^2)]) (((q (μ^2) ℧ cr[τ] (Sqrt[-(℧^2)+2 cr[τ]]) (vT^2))/(((a^2)+(℧^2)+(-2+cr[τ]) cr[τ]) (Sqrt[(a^2)+((cr[τ])^2)])))+((((a^2)+(a^2) Cos[2 cθ[τ]]+2 ((cr[τ])^2)) cr'[τ])/(2 ((a^2)+((cr[τ])^2))))+(((Sqrt[-(℧^2)+2 cr[τ]]) (ct'[τ]-a ((Sin[cθ[τ]])^2) (cφ)'[τ]))/(Sqrt[(a^2)+((cr[τ])^2)]))))/(Sqrt[(((a^2)+(a^2) Cos[2 cθ[τ]]+2 ((cr[τ])^2))/((a^2)+((cr[τ])^2)))]))&&
vΘ==(Sqrt[(a^2) ((Cos[cθ[τ]])^2)+((cr[τ])^2)]) (Sqrt[1-(μ^2) (vT^2)]) (cθ)'[τ]&&
vΦ==(((Sin[cθ[τ]])^2) (Sqrt[1-(μ^2) (vT^2)]) (a q (μ^2) ℧ cr[τ] (Sqrt[(a^2)+((cr[τ])^2)]) (vT^2)-(1/2) a (Sqrt[-(℧^2)+2 cr[τ]]) ((a^2)+(a^2) Cos[2 cθ[τ]]+2 ((cr[τ])^2)) cr'[τ]+(Sqrt[(a^2)+((cr[τ])^2)]) (a ((℧^2)-2 cr[τ]) ct'[τ]+(((((a^2)+((cr[τ])^2)))^2)-(a^2) ((a^2)+(℧^2)-2 cr[τ]+((cr[τ])^2)) ((Sin[cθ[τ]])^2)) (cφ)'[τ])))/((Sqrt[(a^2)+((cr[τ])^2)]) ((a^2) ((Cos[cθ[τ]])^2)+((cr[τ])^2)) (Sqrt[((((Sin[cθ[τ]])^2) (((((a^2)+((cr[τ])^2)))^2)-(a^2) ((a^2)+(℧^2)-2 cr[τ]+((cr[τ])^2)) ((Sin[cθ[τ]])^2)))/((a^2) ((Cos[cθ[τ]])^2)+((cr[τ])^2)))]))&&
vT==(Sqrt[(vR^2)+((vΘ)^2)+((vΦ)^2)]), {vR, vΘ, vΦ , vT}, Reals]; (* Lokalgeschwinditkeitssolver *)

vr[τ_]:=vR/.vf[τ][[1]]; (* Radialgeschwindigkeit *)
vθ[τ_]:=vΘ/.vf[τ][[1]]; (* Longitudinalgeschwindigkeit *)
vφ[τ_]:=vΦ/.vf[τ][[1]]; (* Latitudinalgeschwindigkeit *)
vt[τ_]:=vT/.vf[τ][[1]]; (* Gesamtgeschwindigkeit *)

Plot[Chop@ct[τ], {τ,0,plunge}, Frame->True, 
ImageSize->420, ImagePadding->{{100,1},{20,10}},
PlotPoints->50, MaxRecursion->2, PlotRange->All]

Plot[Chop@cr[τ], {τ,0,plunge}, Frame->True, 
ImageSize->420, ImagePadding->{{100,1},{20,10}},
PlotPoints->50, MaxRecursion->2, PlotRange->All]

Plot[180/π Chop@cθ[τ], {τ,0,plunge}, Frame->True, 
ImageSize->420, ImagePadding->{{100,1},{20,10}},
PlotPoints->50, MaxRecursion->2, PlotRange->All]

Plot[180/π Chop@cφ[τ], {τ,0,plunge}, Frame->True, 
ImageSize->420, ImagePadding->{{100,1},{20,10}},
PlotPoints->50, MaxRecursion->2, PlotRange->All]

Plot[Chop@vr[τ], {τ,0,plunge}, Frame->True, 
ImageSize->420, ImagePadding->{{100,1},{20,10}},
PlotPoints->50, MaxRecursion->2, PlotRange->All]

Plot[Chop@vθ[τ], {τ,0,plunge}, Frame->True, 
ImageSize->420, ImagePadding->{{100,1},{20,10}},
PlotPoints->50, MaxRecursion->2, PlotRange->All]

Plot[Chop@vφ[τ], {τ,0,plunge}, Frame->True, 
ImageSize->420, ImagePadding->{{100,1},{20,10}},
PlotPoints->50, MaxRecursion->2, PlotRange->All]

Plot[Chop@vt[τ], {τ,0,plunge}, Frame->True, 
ImageSize->420, ImagePadding->{{100,1},{20,10}},
PlotPoints->50, MaxRecursion->2, PlotRange->All]

plunge

Table[Chop@Evaluate[r[τ]/.sol][[1]],{τ,0,plunge,0.02}]

Quit[]