I have a question about the difference in the solution between `DSolve`

and `NDSolve`

.

I want to solve the Friedmann equation of

\begin{align} 3\left(\frac{\dot{a}^2}{a^2}+\frac{k}{a^2}\right)=8\pi G\rho(a)\,, \end{align}

where $ a(t)$ is the scale factor, $ k=1,0,-1$ , $ G$ is Newton’s gravitational constant, and $ \rho$ is the energy density. For a Universe composed of radiation, I know that $ \rho\sim a^{-4}$ . For $ k=1$ (and setting $ G=1$ ), I can solve for $ a(t)$ with Mathematica:

` DSolve[{3 (D[a[t], t]^2/a[t]^2 + 1/a[t]^2) == 8 Pi/a[t]^4, a[0] == 1}, a[t], t] `

giving me a solution of:

`{{a[t] -> Sqrt[3 - 2 Sqrt[3 (-3 + 8 \[Pi])] t - 3 t^2]/Sqrt[3]}, {a[t] -> Sqrt[3 + 2 Sqrt[3 (-3 + 8 \[Pi])] t - 3 t^2]/Sqrt[3]}} `

Plotting the second solution, I get:

`Plot[Sqrt[3 + 2 Sqrt[3 (-3 + 8 \[Pi])] t - 3 t^2]/Sqrt[2], {t, 0, 6}] `

which is exactly what I want.

Now, solving the same equation with `NDSolve`

:

`sol = NDSolve[{3 (D[a[t], t]^2/a[t]^2 + 1/a[t]^2) == 8 Pi/a[t]^4, a[0] == 1}, a, {t, 0, 6}] `

I get the following errors:

`NDSolve::ndsz: At t == 0.17823474500439837, step size is effectively zero; singularity or stiff system suspected.`

`NDSolve::mxst: Maximum number of 799782 steps reached at the point t == 5.298557045191269.`

Plotting the solution, I get

`Plot[Evaluate[a[t] /. sol[[2]]], {t, 0, 5.3}] `

We see that the first half of the plot looks similar to the analytical solution, but then nothing is plotted. I decided to plot the real part of the solution and I found

while the imaginary part of the solution gave this plot:

Why are there differences between `DSolve`

and `NDSolve`

? I tried adding the solution

` Method -> {"MethodOfLines", "DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 1}} `

as suggested here, but it didn’t fix the problem. I’ve also tried increasing `MaxSteps`

to more than 100000 as well as changing the `WorkingPrecision`

, `PrecisionGoal`

, and `AccuracyGoal`

. How do I make my numerical solution agree with my analytical solution?