I have the following program to solve for the neutron to baryon ratio during BBN. I've replaced many of the derived values with constants in order to focus on the problem:

```
ClearAll[rate, hubble, equation, initialConditions, fraction]
rate[x_] := (0.28*(12 + 6*x + x^2))/x^5
hubble[x_] := 1.73*^-13/Sqrt[x]
equation := Derivative[1][fraction][x] == (rate[x]/(x*hubble[x]))*((1 - fraction[x])/E^x - fraction[x]); 
initialConditions := fraction[0.12] == 0.46
fraction[x_] = fraction[x] /. NDSolve[{equation, initialConditions}, fraction[x], {x, 0.3, 100}]
LogLogPlot[2*fraction[x], {x, 0.65, 65}]

```

In the domain where $x\approx1$, this ODE goes berserk with these error messages:

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

When I try to plot the ODE solution, I get:

`$IterationLimit::itlim: Iteration limit of 4096 exceeded.`

Here is an image of the progress of the ODE solution as a function of **hubble**. I would expect the solution to the ODE to have the shape of the fifth plot in this series as **hubble** approaches zero.

[![enter image description here](https://i.stack.imgur.com/Ygso9.png)](https://i.stack.imgur.com/Ygso9.png)

I've no idea what this means or how to fix it. It appears to be an issue with precision. How can I fix this? My target is to get this section of code working when $hubble=10^{-13}$.