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

Clear[ρ0, rs]

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

outgoingxderout[x_, τ_]:= 1-((2 (Sqrt[((2 π)/3)]) (rs^(3/2)) (Sqrt[ρ0]))/(Sqrt[x]))
outgoingxderin[x_, τ_] :=((-3+2 (Sqrt[6 π]) x (Sqrt[ρ0])+3 (Sqrt[6 π]) (Sqrt[ρ0]) τ)/(-3+3 (Sqrt[6 π]) (Sqrt[ρ0]) τ))
ingoingxderout[x_, τ_] := -1-((2 (Sqrt[((2 π)/3)]) (rs^(3/2)) (Sqrt[ρ0]))/(Sqrt[x]))
ingoingxderin[x_, τ_] := ((3+2 (Sqrt[6 π]) x (Sqrt[ρ0])-3 (Sqrt[6 π]) (Sqrt[ρ0]) τ)/(-3+3 (Sqrt[6 π]) (Sqrt[ρ0]) τ))
outgoingxder[x_, τ_] := If[outside[x, τ], outgoingxderout[x, τ], outgoingxderin[x, τ]] 
ingoingxder[x_, τ_] := If[outside[x, τ], ingoingxderout[x, τ], ingoingxderin[x, τ]] 
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, τ], xderout[x, τ], xderin[x, τ]]

(*Verification that the piecewise solution matches the solution where the ingoing ray is simply taken to travel from -x *)
solin = NDSolve[{x'[τ] == ingoingxder[x[τ], τ], x[-0.505] == 0}, x[τ], {τ, -0.505, -10.0}]
solout = NDSolve[{x'[τ] == outgoingxder[x[τ], τ], x[0.0] == 2 m + .01}, x[τ], {τ, -10.0, 10.0}]
solout2 = NDSolve[{x'[τ] == outgoingxder[x[τ], τ], x[0.0] == 2 m + .1}, x[τ], {τ, -10, 10.0}]
solin2 = NDSolve[{x'[τ] == ingoingxder[x[τ], τ], x[-0.627] == 0}, x[τ], {τ, -0.627, -10.0}]
solout3 = NDSolve[{x'[τ] == outgoingxder[x[τ], τ], x[0.0] == 2 m+.001}, x[τ], {τ, -10, 10.0}]
solin3 = NDSolve[{x'[τ] == ingoingxder[x[τ], τ], x[-0.504] == 0}, x[τ], {τ, -0.504, -10.0}]
sol4 = NDSolve[{x'[τ] == xder[x[τ], τ], x[0.0] == 2 m+.001}, x[τ], {τ, -10, 10.0}]
Plot[{If[τ < -0.505, x[τ] /. solin[[1]], x[τ] /.solout[[1]]], If[τ < -0.627, x[τ] /. solin2[[1]], x[τ] /.solout2[[1]]], If[τ < -0.504, x[τ] /. solin3[[1]], x[τ] /.solout3[[1]]], Abs[x[τ]] /. sol4},  {τ, -10.0, 10.0}, PlotRange -> Full]

(*The fourth graph is on top of the others, so we don't need to consider the piecewise solutions any more *)

ρ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}, PlotRange -> Full]

{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]]}

Plot[{ Abs[x[τ]] /. sol, Abs[x[τ]] /. sol1, Abs[x[τ]] /. sol3, Abs[x[τ]] /. sol4, xout[rs, τ]}, {τ, -2, 2}, PlotRange -> Full]