Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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!

share|improve this question
Hey, welcome! Try to find a minimal example that illustrate your problem, so that it shows your effort, is less work for the people trying to help you, and is more general and helps future visitors more – Rojo Feb 22 '13 at 18:17
Hiya, sorry about that. I've condensed the problem down a bit...thanks for pointing that out! – Gokotai Feb 22 '13 at 18:28
You need to provide numerical values for all the constants too. – b.gatessucks Feb 22 '13 at 18:54
done! sorry about that too... – Gokotai Feb 22 '13 at 18:55
up vote 3 down vote accepted

You could do:

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

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

Mathematica graphics

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

share|improve this answer
as a side note, you could also use NDSolveValue that would avoid that ... /. ScaleFactor[[1]] business. – user21 Feb 23 '13 at 8:45
@ruebenko Not me :) ... still on v.8 – Dr. belisarius Feb 23 '13 at 8:48
ah, OK, then as an idea for the OP, and maybe a small incentive to upgrade ;-) – user21 Feb 23 '13 at 8:52

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.