The following is my code:
a := 3.24077*10^-20 (* km \[Rule] Mpc *)
b := 3.16888*10^-14 (* s \[Rule] MYear *)
c := a/b*(2.99792*10^5) (*Mpc/ MYear*)
H0 := a/b*71 (*1/MYear*)
G := a^3/b^2*6.67398*10^-20 (*Mpc^3/(Kg*MYear^2)*)
\[Rho]crit1 := 3/(8 \[Pi]*G)*(H0)^2
M1[r_] := (4 \[Pi] )/3*\[Rho]crit1*(r)^3
E1[r_] := 0
ScaleFactor =
NDSolve[{Sqrt[R[r, t]] (D[R[r, t], t]) ==
Sqrt[(2*G*M1[r])/c + 2 c*E1[r]*R[r, t]], R[r, 0] == r/1000},
R, {r, 1, 10000}, {t, 1, 10000}] // FullSimplify;
Plot3D[Evaluate[R[r, t] /. %], {r, 1, 10000}, {t, 1, 10000}]
ParticleHorizon[r_, g_] :=
NIntegrate[(c Sqrt[1 + 2 E1[r]])/
D[Evaluate[R[r1, t] /. ScaleFactor, r1 -> r], r], {t, 0, g}]
Plot3D[ParticleHorizon[r,g], {r, 1, 10}, {g, 1, 10}]
My problem is that once I've solved the differential equation to get R[r,t]
(which it solves without a problem, giving me an interpolating function), I have trouble differentiating the result with respect to r
, to be used in my ParticleHorizon
function. I think it might be because the differential operator uses DSolve
or NDSolve
intrinsically.
Here, a,b,c,f,G
and lambda
are constants and M1[r]
is a boundary condition that we need to solve the differential equation for R[r,t]
. E1[r]
has been set to 0 for now for simplicity.
Any ideas on how to proceed?
Also, the rule to change r1->r
was just me trying to see if it works, which it doesn't...
Thanks in advance!