rout[x_, τ_] := (x^(3/2) + Sqrt[6 π ρ0] τ rs^(3/2))^(2/3)
rin[x_, τ_] := x/(1 - Sqrt[6 π ρ0] τ)^(2/3)
xin[r_, τ_] := r (1 - Sqrt[6 π ρ0] τ)^(2/3)
xout[r_ , τ_] := (r^(3/2) - Sqrt[6 π ρ0] τ rs^(3/2))^(2/3)
outside[x_, τ_] := (Abs[x]^(3/2) > rs^(3/2) (1 - Sqrt[6 π ρ0] τ))

(*Metric outside in terms of x *)
metricout = -dτ^2 + (D[rout[x, τ], x] dx  + D[rout[x, τ], τ] dτ)^2/(1 - Sqrt[6 π ρ0] τ rs^(3/2)/rout[x, τ]^(3/2))^(2/3) //Simplify//Expand

Coefficient[metricout, -dτ^2] //Simplify//PowerExpand
Coefficient[metricout, dx dτ]/2 //PowerExpand
Coefficient[metricout, dx^2] //Simplify //PowerExpand

gttout[x_, τ_] :=1-((8 π (rs^3) ρ0)/(3 x))
gxtout[x_, τ_] := ((2 (Sqrt[((2 π)/3)]) (rs^(3/2)) (Sqrt[ρ0]))/(Sqrt[x]))
gxxout[x_, τ_] := 1

(* Metric Inside in terms of x*)
metricin= -dτ^2 + (D[rin[x, τ], x] dx  + D[rin[x, τ], τ] dτ)^2(1 - Sqrt[6 π ρ0] τ)^(4/3) //Simplify//Expand

Coefficient[metricin, -dτ^2] //Simplify//PowerExpand
Coefficient[metricin, dx dτ]/2 //PowerExpand
Coefficient[metricin, dx^2]//PowerExpand

gttin[x_, τ_] := ((3-8 π (x^2) ρ0-6 (Sqrt[6 π]) (Sqrt[ρ0]) τ+18 π ρ0 (τ^2))/(3 (((-1+(Sqrt[6 π]) (Sqrt[ρ0]) τ))^2)))
gxtin[x_, τ_] := ((2 (Sqrt[((2 π)/3)]) x (Sqrt[ρ0]))/(1-(Sqrt[6 π]) (Sqrt[ρ0]) τ))
gxxin[x_, τ_] := 1

oxout=(-gxtout[x, τ] + Sqrt[gxtout[x, τ]^2 +  gttout[x, τ] gxxout[x, τ]])/gxxout[x, τ]//Simplify//PowerExpand//Simplify
 oxin = (-gxtin[x, τ] + Sqrt[gxtin[x, τ]^2 +  gttin[x, τ] gxxin[x, τ]])/gxxin[x, τ] //Simplify//PowerExpand//Simplify
 ixout = (-gxtout[x, τ] - Sqrt[gxtout[x, τ]^2 +  gttout[x, τ] gxxout[x, τ]])/gxxout[x, τ] //Simplify//PowerExpand//Simplify
 ixin = (-gxtin[x, τ] - Sqrt[gxtin[x, τ]^2 +  gttin[x, τ] gxxin[x, τ]])/gxxin[x, τ] //Simplify//PowerExpand//Simplify

xderout[x_, τ_] := 1-((2 (Sqrt[((2 π)/3)]) (rs^(3/2)) (Sqrt[ρ0]))/(Sign[x](Sqrt[Abs[x]]))) (*This tries to capture the change in sign of the square root via the change in sign in x*)
xderin[x_, τ_] := ((-3+2 (Sqrt[6 π]) x (Sqrt[ρ0])+3 (Sqrt[6 π]) (Sqrt[ρ0]) τ)/(-3+3 (Sqrt[6 π]) (Sqrt[ρ0]) τ)) (*This is also a combined function that captures the change in sign of the square root via the change in sign in x*)
xder[x_, τ_] :=  If[outside[x, τ],
	If[Abs[x] < .001, 0,xderout[x, τ]], 
	xderin[x, τ]] 
(*This sets a cutoff so that if the ray is outside and also x is very small (which means that it is already trapped inside and heading to the singularity, we set the derivative to zero to avoid hitting the singularity and having the equation blow up.*)

ρ0 = (1/ (6 π));
rs = 1;
m =(4/3) π rs^3 ρ0
sol= NDSolve[{x'[τ] == xder[x[τ], τ], x[1.0] == 2 m+.001}, x, {τ, -30.0, 100.0}]
sol1= NDSolve[{x'[τ] == xder[x[τ], τ], x[1.0] == 2 m+.01}, x, {τ, -30, 100.0}]
sol3= NDSolve[{x'[τ] == xder[x[τ], τ], x[1.0] == 2 m}, x, {τ, -30.0, 100.0}]
sol4= NDSolve[{x'[τ] == xder[x[τ], τ], x[1.0] == 2 m+.00001}, x, {τ, -30, 100.0}]

Plot[{ Abs[x[τ]] /. sol, Abs[x[τ]] /. sol1, Abs[x[τ]] /. sol3, Abs[x[τ]] /. sol4, xout[rs, τ]}, {τ, -30, 30.0}, PlotRange -> Full]
(*Both sol3 and sol4 are inside the shell at the beginning, whereas sol1 and sol2 are outside. Note how tiny tiny differences at one point of time, get blown up in large differences at another point of time *)

{x[-30]/.sol[[1]], x[-30]/.sol1[[1]], x[-30]/.sol3[[1]], x[-30]/.sol4[[1]]}
{x[30]/.sol[[1]], x[30]/.sol1[[1]], x[30]/.sol3[[1]], x[30]/.sol4[[1]]}