ClearAll[aS,aK,ro,rmax ]
Needs["PlotLegends`"]
aS=0;
aK=1;
(*rK=1+(Sqrt[1-(aK^2)]);*)
ro=0;
rmax=25;
Sig[r_,a_]:=(r^2)+(a^2);
T[r_,a_,l_]:=(r^2)+(a^2)-l a;
delta[r_,a_]:=(r^2)+(a^2)-2r;
grr[r_,a_]:=((Sig[r,a])/(delta[r,a]));
gtt[r_,a_]:=-(1-((2r)/(Sig[r,a])));
Vr[r_,a_,l_]:=(((T[r,a,l]))^2)-delta[r,a] ((r^2)+(((l-a))^2));
drdtau[r_,a_,l_]:=If[l>0,((Sqrt[Vr[r,a,l]])/(Sig[r,a])),-((Sqrt[Vr[r,a,l]])/(Sig[r,a]))];
drdtaulp[r_,a_,l_]:=((Sqrt[Vr[r,a,l]])/(Sig[r,a]));
drdtaulm[r_,a_,l_]:=-((Sqrt[Vr[r,a,l]])/(Sig[r,a]));
dtdtau[r_,a_,l_]:=(1/(Sig[r,a]))(-a (a-l)+((r^2)+(a^2))((T[r,a,l])/(delta[r,a])));
lapse[r_,a_,l_]:=(1/(dtdtau[r,a,l]));
drdt[r_,a_,l_]:=((drdtau[r,a,l])/(dtdtau[r,a,l]));
(*the following should be compared to the left panels of Figs. 2 and 3 of http://journals.aps.org/prl/pdf/10.1103/PhysRevLett.103.111102, which in turn should be compared to good old-fashioned physical intuition*)
Plot[{drdt[r,aS,2],drdt[r,aS,4],drdt[r,aS,5]},{r,ro,rmax},PlotRange->{{ro,rmax},{-0.4,0.4}},PlotStyle->{{Red,Dashed,Thick},{Blue,Thick},{Black,Dashing[{.02,.02}],Thick}},Frame->True,FrameLabel->{"r[M]","dr/dt"},LabelStyle->Directive[Black,FontSize->24], PlotLegend->{Style["l=2, Schw",24],Style["l=4, Schw",24],Style["l=5, Schw",24]},LegendPosition->{.2,-.2},LegendSize->{.8,.4},LegendShadow->None,AxesOrigin->{0,0}]
Plot[{drdt[r,aK,0.4],drdt[r,aK,2],drdt[r,aK,3]},{r,ro,rmax},PlotRange->{{ro,rmax},{-0.4,0.4}},PlotStyle->{{Red,Dashed,Thick},{Blue,Thick},{Black,Dashing[{.02,.02}],Thick}},Frame->True,FrameLabel->{"r[M]","dr/dt"},LabelStyle->Directive[Black,FontSize->24], PlotLegend->{Style["l=0.4, Kerr",24],Style["l=2, Kerr",24],Style["l=3, Kerr",24]},LegendPosition->{.2,-.2},LegendSize->{.8,.4},LegendShadow->None,AxesOrigin->{0,0}]