Problems with doing calculus on Interpolating functions

Question

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!

Answer 1

You could do:

ParticleHorizon[r_, g_] := 
   NIntegrate[(c Sqrt[1 + 2 E1[r]])/Evaluate[D[R[r, t] /. ScaleFactor[[1]]], r], {t, 0, g}]

ListPlot3D[
        Flatten[Table[{r, g, ParticleHorizon[r, g]}, {r, 1, 10}, {g, 1, 10}], 1], 
        ColorFunction -> "SouthwestColors", Mesh -> True, 
        MeshFunctions -> {#3 &}, 
        PlotStyle -> Directive[Opacity[0.8], Specularity[White, 50]]]

I used ListPlot3D just because it is faster than Plot3D[]