# Numerical artifact in solution of system of ODEs

**URL:** <https://discourse.julialang.org/t/numerical-artifact-in-solution-of-system-of-odes/130941>\
**Category:** Numerics\
**Created:** [July 22, 2025, 11:47am UTC](https://discourse.julialang.org/t/numerical-artifact-in-solution-of-system-of-odes/130941 "2025-07-22T11:47:29Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![Gako](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gako/32/217772_2.png) [@Gako](https://discourse.julialang.org/u/Gako)\
**Post date:** [July 22, 2025, 11:47am UTC](https://discourse.julialang.org/t/numerical-artifact-in-solution-of-system-of-odes/130941/1 "2025-07-22T11:47:29Z")

</div>

I am currently solving the following system of ODEs:

\begin{cases} \ddot\phi=-3H\dot{\phi}-\phi\\ \dot{a}=Ha\\ \ddot\sigma=-H\dot{\sigma}-\frac{m\_{Eff}^2}{a^2}\sigma\\ \ddot{\rho}\_0=\frac{4\rho\_2}{a^2}-H\dot{\rho}\_{0}-\frac{4\omega^2\rho\_{0}}{a^2} \end{cases}

where the last equation is solved for different k\in[1, 10^4] and

\begin{cases} H=\sqrt{\frac{\dot{\phi }^2+\phi ^2}{6}}\\ R=\dot{\phi }^2-2\phi^2\\ m\_{Eff}^2=-a^2R\left(\xi-\frac{1}{6} \right)\\ \omega ^2=k^2+m\_{Eff}^2\\ \rho \_2=\frac{1}{2\rho \_{0}}\left(0.25 + \omega^2\rho \_0^2+\frac{\dot{\rho }\_{0}^2}{4}\right).\end{cases}

Here’s an MWE of it (with \xi = \frac{1}{5}):

```julia
using DifferentialEquations
using Plots

function solve_equations()
    kavlength = 599
    k = Vector(logrange(1, 10^4, kavlength + 1))
    kav = (k[1:end - 1] + k[2:end]) / 2
    iconds = initialConditions(kav, kavlength)
    tspan = (0.0, 1.5)
    teval = range(tspan[1], tspan[2], 501)
    ddrho0 = zeros(kavlength)
    omega2 = zeros(kavlength)

    prob = ODEProblem(diffEquations!, iconds, tspan, (kav, kavlength, ddrho0, omega2), saveat=teval)
    sol = solve(prob, Vern9(), abstol=1e-15, reltol=1e-15)
    # savefig(heatmap(sol.t, kav, log10.(abs.(sol[6:5 + kavlength, :])), yscale=:log10, dpi=600, title="\\delta\\rho_0", xlabel="t", ylabel="k"), "heatmap.svg")
end

function initialConditions(kav, kavlength)
    rho0 = ones(kavlength) ./ (2 .* sqrt.(16 .+ kav.^2))
    drho0 = zeros(size(rho0))
    arr = [15, -0.8, 1, 0.001, 0]
    initial = vcat(arr, rho0, drho0)
    return initial
end

function ddrho0!(a, H, omega2, rho0, drho0, ddrho0)
    for i in eachindex(omega2)
        rho2 = (0.25 + omega2[i] * rho0[i]^2 + drho0[i]^2 / 4) / (2 * rho0[i])
        ddrho0[i] = 4 * rho2 / a^2 - H * drho0[i] - 4 * omega2[i] * rho0[i] / a^2
    end
    nothing
end

function diffEquations!(du, u, args, t)
    kav, kavlength, ddrho0, omega2 = args

    infl, dinfl, a, sigma, dsigma = @view u[1:5]
    rho0 = @view u[6:6 + kavlength - 1]
    drho0 = @view u[6 + kavlength:6 + 2 * kavlength - 1]

    H = sqrt((dinfl^2 + infl^2) / 6)
    da = a * H
    R = dinfl^2 - 2 * infl^2
    
    m2Eff = -a^2 * R / 30
    @. omega2 = kav^2 + m2Eff
    ddrho0!(a, H, omega2, rho0, drho0, ddrho0)
    ddinfl = -3 * H * dinfl - infl
    ddsigma = - H * dsigma - m2Eff * sigma / a^2

    du[1] = dinfl
    du[2] = ddinfl
    du[3] = da
    du[4] = dsigma
    du[5] = ddsigma
    du[6:5 + kavlength] .= drho0
    du[6 + kavlength:5 + 2 * kavlength] .= ddrho0
    nothing
end

solve_equations()

```

The solver actually errors at t\approx1,35 with the error message

```julia
Warning: At t=1.3482231005865255, dt was forced below floating point epsilon 2.220446049250313e-16, and step error estimate = 14.835297241231595. Aborting. There is either an error in your model specification or the true solution is unstable (or the true solution can not be represented in the precision of Float64).

```

From the solution one can create a heatmap of log\_{10}(\delta\rho\_0) which looks like this:

 ![Heatmap of log10(deltarho0)](https://global.discourse-cdn.com/julialang/original/3X/d/f/dfa6985e96a3dc691b9784a2ef24b2ff049b5c89.png)  
There a numerical artifact can be seen starting from around t=0,6. I believe another one is appearing from the point the solver fails as well. The solver already had issues getting over the first line, slowing down considerably. Changing the density of k values doesn’t change the line and changing k\_{min} and k\_{max} doesn’t change its position either, just the part that can be seen. The solver tolerances are already at 1e-15, changing them one way or another doesn’t make a difference either. Digging into the values of the solution the only notable thing is a seemingly large gap in the derivative of \delta\rho\_0 at where the line is. Here’s an example of the derivative of \delta\rho\_0 evaluated for k\approx1:

| t | derivative of deltarho0 |
| --- | --- |
| 0.564 | -0.140 |
| 0.567 | -0.137 |
| 0.57 | -0.134 |
| 0.573 | -0.130 |
| 0.576 | -0.127 |
| 0.579 | -0.123 |
| 0.582 | -0.116 |
| 0.585 | 0.114 |
| 0.588 | 0.116 |
| 0.591 | 0.116 |
| 0.594 | 0.115 |
| 0.597 | 0.114 |
| 0.6 | 0.112 |
| 0.603 | 0.110 |

How would a jump (and line) like this appear?

---

<div class="post-metadata">

**Author:** ![John\_Gibson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/john_gibson/32/5321_2.png) [@John\_Gibson](https://discourse.julialang.org/u/John_Gibson)\
**Post date:** [July 22, 2025, 12:11pm UTC](https://discourse.julialang.org/t/numerical-artifact-in-solution-of-system-of-odes/130941/2 "2025-07-22T12:11:32Z")

</div>

What are k and \xi in relation to the differential equation? They do not appear in the ODE.

---

<div class="post-metadata">

**Author:** ![Gako](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gako/32/217772_2.png) [@Gako](https://discourse.julialang.org/u/Gako)\
**Post date:** [July 22, 2025, 12:30pm UTC](https://discourse.julialang.org/t/numerical-artifact-in-solution-of-system-of-odes/130941/3 "2025-07-22T12:30:57Z")

</div>

They are effectively constants. \xi is a parameter of the system, in the MWE I’ve set it to 0,2. k is also a constant but the fourth DE is evaluated with different values of k within the range k\in [1, 10^4].

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [July 22, 2025, 12:32pm UTC](https://discourse.julialang.org/t/numerical-artifact-in-solution-of-system-of-odes/130941/4 "2025-07-22T12:32:52Z")

</div>

In the code, \xi = 1/5 = 0.2 (using the English decimal point…), thus, \xi - 1/6 = 1/30.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [July 22, 2025, 12:57pm UTC](https://discourse.julialang.org/t/numerical-artifact-in-solution-of-system-of-odes/130941/5 "2025-07-22T12:57:08Z")

</div>

Are you sure it’s not a singularity or instability in your equation? e.g. \rho\_2 blows up whenever \rho\_0 passes through zero, which could account for the apparent discontinuity in \dot{\rho}\_0.

---

<div class="post-metadata">

**Author:** ![Gako](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gako/32/217772_2.png) [@Gako](https://discourse.julialang.org/u/Gako)\
**Post date:** [July 22, 2025, 1:10pm UTC](https://discourse.julialang.org/t/numerical-artifact-in-solution-of-system-of-odes/130941/6 "2025-07-22T13:10:14Z")

</div>

I’m quite confident in that (albeit less so after seeing it staying at the same place with different k\_{min} and k\_{max}). \rho\_0 shouldn’t, and indeed doesn’t, go negative, at least on the points that were saved. It is possible, however, that I’ve misstyped something in the equations that’s creating this effect.

---

<div class="post-metadata">

**Author:** ![AntonReinhard](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antonreinhard/32/211306_2.png) [@AntonReinhard](https://discourse.julialang.org/u/AntonReinhard)\
**Post date:** [July 22, 2025, 5:24pm UTC](https://discourse.julialang.org/t/numerical-artifact-in-solution-of-system-of-odes/130941/7 "2025-07-22T17:24:07Z")

</div>

I’m not 100% sure if this is a viable idea because I don’t know much about solving systems of ODEs, but you could try using my new package PrecisionCarriers.jl to see if your functions maybe become numerically unstable around that line. There are a few subtractions in your functions that could lead to catastrophic cancellations for example.

I already tried using your code myself and solving the entire thing using the package, but it seems that DifferentialEquations.jl casts the values to Float64 somewhere (in the Vern9 cache?), which unfortunately gets rid of the carried precision. I’m not sure if this is something I could fix for PrecisionCarriers.

---

<div class="post-metadata">

**Author:** ![Sevi](https://avatars.discourse-cdn.com/v4/letter/s/c67d28/32.png) [@Sevi](https://discourse.julialang.org/u/Sevi)\
**Post date:** [July 22, 2025, 7:09pm UTC](https://discourse.julialang.org/t/numerical-artifact-in-solution-of-system-of-odes/130941/8 "2025-07-22T19:09:48Z")

</div>

Could you break the code down to a smaller example? If I understand you correctly, k is a fixed parameter in the set of ODEs you solve, but in the code it looks like you solve the ODEs for many different values of k simultaneously in one big equation.

If these equations all decouple, I’d wager it’s faster and more stable to solve them individually with one problem per value of k (perhaps with an `EnsembleProblem`). Otherwise, you are constraining the solver for all values of k to the equations with the “most unfavorable” k-values.

In any case, a simpler MWE would make it easier to dig further into the issue. Unless you are actually solving a discretized PDE, in which case other things might be going on.

---

<div class="post-metadata">

**Author:** ![apo383](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/apo383/32/11272_2.png) [@apo383](https://discourse.julialang.org/u/apo383)\
**Post date:** [July 22, 2025, 7:38pm UTC](https://discourse.julialang.org/t/numerical-artifact-in-solution-of-system-of-odes/130941/9 "2025-07-22T19:38:26Z")

</div>

What’s going on with a? It looks unstable, because \dot{a} = Ha and H is a positive square root.

All sorts of weird things happen, for example m^2\_{Eff}=-a^2... should get big (not sure if positive or negative), and so \omega^2 should get big, and so \rho\_2 should get big even if \rho\_0 doesn’t pass through zero?

Several differential right-hand-sides multiply and then divide by a^2 in code, which doesn’t seem like good policy as a blows up. All this could be wrong depending on the constants, but it’s hard to see what’s going on, hence questions.

The dynamics look suspicious to me. I wonder if there is a simplification of the equations that eliminates need to multiply and divide by an unstable quantity.

As with @Sevi, I recommend ditching the sweep of k, since you get the artifact for all small values seen in the plot (and probably larger ones as well). Just choose one k, and plot all the states over time to see what they do. And straighten out the situation with a.

---

<div class="post-metadata">

**Author:** ![Gako](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gako/32/217772_2.png) [@Gako](https://discourse.julialang.org/u/Gako)\
**Post date:** [July 23, 2025, 6:44am UTC](https://discourse.julialang.org/t/numerical-artifact-in-solution-of-system-of-odes/130941/10 "2025-07-23T06:44:10Z")

</div>

Thank you for your replies. Unfortunately, as so it often happens, the problem was between the keyboard and the chair. While writing the MWE for a single k I went over the equations once more and noticed that I had forgotten a single a^2 term from the equation for \rho\_2 (which should’ve been \rho\_2 = \frac{1}{2\rho\_0}\left(0.25+\omega^2\rho\_0^2+\frac{a^2\dot\rho\_0^2}{4}\right)). With this correction the “artifact” disappeared and @stevengj was indeed correct, it was an instability of the system.

Notwithstanding, I edited the calculations slightly based on @apo383 's suggestions, mainly changing m\_{Eff}^2 for \frac{m\_{Eff}^2}{a^2} = -R\left(\xi - \frac{1}{6}\right) and \omega^2 for \frac{\omega^2}{a^2} = \frac{k^2}{a^2} + \frac{m\_{Eff}^2}{a^2} and merging the equations for \ddot\rho\_0 and \rho\_2 to \ddot\rho\_0 = \frac{1}{2\rho\_0}\left(\dot\rho\_0^2 + \frac{1}{a^2}\right) - H\dot\rho\_0 -2\rho\_0\frac{\omega^2}{a^2}. Hopefully they’ll help with the stability with larges values of t.

 ![rho0 without a huge rift in the middle](https://global.discourse-cdn.com/julialang/original/3X/3/8/388a0b90074bc71f54cd6bb5cb40d0b7bc4008b2.png)

---

<div class="post-metadata">

**Author:** ![apo383](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/apo383/32/11272_2.png) [@apo383](https://discourse.julialang.org/u/apo383)\
**Post date:** [July 23, 2025, 8:39pm UTC](https://discourse.julialang.org/t/numerical-artifact-in-solution-of-system-of-odes/130941/11 "2025-07-23T20:39:36Z")

</div>

It looks like a good idea to get rid of \rho\_2. An additional suggestion is to replace dividing by a^2 with multiplying by b^2, where \dot{b}=-Hb. The problem with a is that it’s unstable and heading to an overflow error. Whereas b is decaying exponentially, and will gracefully underflow to zero.

Perhaps \rho\_0 is unstable so you’re headed for trouble no matter what, but still seems best to isolate a single instability, otherwise it’s a race to see which one overflows first.

In general I also find it helpful to set expectations with a toy problem. The fact that you complained about the line artifact suggests you didn’t have an expectation on the overall behavior. A toy problem is one with trivial parameter values, where you set some to zero or infinity or something easy to reason with. Right now it’s hard to reason about your system because it’s high order and coupled. Perhaps \phi and \sigma are supposed to be stable, so can you set them to zero or constant? If you thus decouple \phi\_0, it may be evident whether it should be unstable.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 24, 2025, 9:46am UTC](https://discourse.julialang.org/t/numerical-artifact-in-solution-of-system-of-odes/130941/12 "2025-07-24T09:46:29Z")

</div>

Did you make u0 and t be the type? Or just u0? Because it is really well tested that DiffEq only casts to types that the user gives, and I have seen people post an issue about this where they only change the ODE state but keep a time type just Float64. Those two types are kept separate because for example u0 can be complex while t cannot. So if you want to do something to trace position, you should put the tracer to both u0 and t and that should work just fine (and is really well tested with things like AD)

---

<div class="post-metadata">

**Author:** ![AntonReinhard](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antonreinhard/32/211306_2.png) [@AntonReinhard](https://discourse.julialang.org/u/AntonReinhard)\
**Post date:** [July 24, 2025, 10:09am UTC](https://discourse.julialang.org/t/numerical-artifact-in-solution-of-system-of-odes/130941/13 "2025-07-24T10:09:46Z")

</div>

Hm, as I said I don’t have experience with DifferentialEquations.jl, but as far as I can tell and debug, every variable is a PrecisionCarrier.

I just took the OP’s code and put precify on everything:

```julia
using DifferentialEquations
using Plots
using PrecisionCarriers

function solve_equations()
    kavlength = 599
    k = precify(Vector(logrange(1, 10^4, kavlength + 1)))
    kav = precify((k[1:end - 1] + k[2:end]) / 2)
    iconds = precify(initialConditions(kav, kavlength))
    tspan = precify((0.0, 1.5))
    teval = range(tspan[1], tspan[2], 501) # already a precified range from tspan
    ddrho0 = precify(zeros(kavlength))
    omega2 = precify(zeros(kavlength))

    prob = ODEProblem(diffEquations!, iconds, tspan, (kav, kavlength, ddrho0, omega2), saveat=teval)
    sol = solve(prob, Vern9(), abstol=precify(1e-15), reltol=precify(1e-15))
    # savefig(heatmap(sol.t, kav, log10.(abs.(sol[6:5 + kavlength, :])), yscale=:log10, dpi=600, title="\\delta\\rho_0", xlabel="t", ylabel="k"), "heatmap.svg")
end

function initialConditions(kav, kavlength)
    rho0 = ones(kavlength) ./ (2 .* sqrt.(16 .+ kav.^2))
    drho0 = zeros(size(rho0))
    arr = [15, -0.8, 1, 0.001, 0]
    initial = vcat(arr, rho0, drho0)
    return initial
end

function ddrho0!(a, H, omega2, rho0, drho0, ddrho0)
    for i in eachindex(omega2)
        rho2 = (0.25 + omega2[i] * rho0[i]^2 + drho0[i]^2 / 4) / (2 * rho0[i])
        ddrho0[i] = 4 * rho2 / a^2 - H * drho0[i] - 4 * omega2[i] * rho0[i] / a^2
    end
    nothing
end

function diffEquations!(du, u, args, t)
    kav, kavlength, ddrho0, omega2 = args

    infl, dinfl, a, sigma, dsigma = @view u[1:5]
    rho0 = @view u[6:6 + kavlength - 1]
    drho0 = @view u[6 + kavlength:6 + 2 * kavlength - 1]

    H = sqrt((dinfl^2 + infl^2) / 6)
    da = a * H
    R = dinfl^2 - 2 * infl^2
    
    m2Eff = -a^2 * R / 30
    @. omega2 = kav^2 + m2Eff
    ddrho0!(a, H, omega2, rho0, drho0, ddrho0)
    ddinfl = -3 * H * dinfl - infl
    ddsigma = - H * dsigma - m2Eff * sigma / a^2

    du[1] = dinfl
    du[2] = ddinfl
    du[3] = da
    du[4] = dsigma
    du[5] = ddsigma
    du[6:5 + kavlength] .= drho0
    du[6 + kavlength:5 + 2 * kavlength] .= ddrho0
    nothing
end

solve_equations()

```

If I do the evil thing and force a `convert(Float64, PrecisionCarrier)` to return a PrecisionCarrier again, I get the following error, because somewhere it seems to cache values into a `Vector{Float64}`:

```julia
TypeError: in typeassert, expected Float64, got a value of type PrecisionCarrier{Float64}

Stacktrace:
  [1] setindex!(A::Vector{Float64}, x::PrecisionCarrier{Float64}, i::Int64)
    @ Base ./array.jl:987
  [2] macro expansion
    @ ~/.julia/packages/DiffEqBase/yhgdI/src/calculate_residuals.jl:107 [inlined]
  [3] macro expansion
    @ ./simdloop.jl:77 [inlined]
  [4] calculate_residuals!
    @ ~/.julia/packages/DiffEqBase/yhgdI/src/calculate_residuals.jl:106 [inlined]
  [5] perform_step!(integrator::OrdinaryDiffEqCore.ODEIntegrator{Vern9{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, true, Vector{PrecisionCarrier{Float64}}, Nothing, PrecisionCarrier{Float64}, Tuple{Vector{PrecisionCarrier{Float64}}, Int64, Vector{PrecisionCarrier{Float64}}, Vector{PrecisionCarrier{Float64}}}, PrecisionCarrier{Float64}, PrecisionCarrier{Float64}, PrecisionCarrier{Float64}, PrecisionCarrier{Float64}, Vector{Vector{PrecisionCarrier{Float64}}}, ODESolution{PrecisionCarrier{Float64}, 2, Vector{Vector{PrecisionCarrier{Float64}}}, Nothing, Nothing, Vector{PrecisionCarrier{Float64}}, Vector{Vector{Vector{PrecisionCarrier{Float64}}}}, Nothing, ODEProblem{Vector{PrecisionCarrier{Float64}}, Tuple{PrecisionCarrier{Float64}, PrecisionCarrier{Float64}}, true, Tuple{Vector{PrecisionCarrier{Float64}}, Int64, Vector{PrecisionCarrier{Float64}}, Vector{PrecisionCarrier{Float64}}}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(diffEquations!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, @Kwargs{saveat::LinRange{PrecisionCarrier{Float64}, Int64}}, SciMLBase.StandardODEProblem}, Vern9{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, OrdinaryDiffEqCore.InterpolationData{ODEFunction{true, SciMLBase.AutoSpecialize, typeof(diffEquations!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, Vector{Vector{PrecisionCarrier{Float64}}}, Vector{PrecisionCarrier{Float64}}, Vector{Vector{Vector{PrecisionCarrier{Float64}}}}, Nothing, OrdinaryDiffEqVerner.Vern9Cache{Vector{PrecisionCarrier{Float64}}, Vector{PrecisionCarrier{Float64}}, Vector{Float64}, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Nothing}, SciMLBase.DEStats, Nothing, Nothing, Nothing, Nothing}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(diffEquations!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, OrdinaryDiffEqVerner.Vern9Cache{Vector{PrecisionCarrier{Float64}}, Vector{PrecisionCarrier{Float64}}, Vector{Float64}, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, OrdinaryDiffEqCore.DEOptions{PrecisionCarrier{Float64}, PrecisionCarrier{Float64}, PrecisionCarrier{Float64}, PrecisionCarrier{Float64}, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{PrecisionCarrier{Float64}, DataStructures.FasterForward}, DataStructures.BinaryHeap{PrecisionCarrier{Float64}, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, LinRange{PrecisionCarrier{Float64}, Int64}, Tuple{}}, Nothing, Float64, Nothing, OrdinaryDiffEqCore.DefaultInit, Nothing}, cache::OrdinaryDiffEqVerner.Vern9Cache{Vector{PrecisionCarrier{Float64}}, Vector{PrecisionCarrier{Float64}}, Vector{Float64}, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, repeat_step::Bool)
    @ OrdinaryDiffEqVerner ~/.julia/packages/OrdinaryDiffEqVerner/rIAEX/src/verner_rk_perform_step.jl:1152
  [6] perform_step!
    @ ~/.julia/packages/OrdinaryDiffEqVerner/rIAEX/src/verner_rk_perform_step.jl:1043 [inlined]
  [7] solve!(integrator::OrdinaryDiffEqCore.ODEIntegrator{Vern9{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, true, Vector{PrecisionCarrier{Float64}}, Nothing, PrecisionCarrier{Float64}, Tuple{Vector{PrecisionCarrier{Float64}}, Int64, Vector{PrecisionCarrier{Float64}}, Vector{PrecisionCarrier{Float64}}}, PrecisionCarrier{Float64}, PrecisionCarrier{Float64}, PrecisionCarrier{Float64}, PrecisionCarrier{Float64}, Vector{Vector{PrecisionCarrier{Float64}}}, ODESolution{PrecisionCarrier{Float64}, 2, Vector{Vector{PrecisionCarrier{Float64}}}, Nothing, Nothing, Vector{PrecisionCarrier{Float64}}, Vector{Vector{Vector{PrecisionCarrier{Float64}}}}, Nothing, ODEProblem{Vector{PrecisionCarrier{Float64}}, Tuple{PrecisionCarrier{Float64}, PrecisionCarrier{Float64}}, true, Tuple{Vector{PrecisionCarrier{Float64}}, Int64, Vector{PrecisionCarrier{Float64}}, Vector{PrecisionCarrier{Float64}}}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(diffEquations!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, @Kwargs{saveat::LinRange{PrecisionCarrier{Float64}, Int64}}, SciMLBase.StandardODEProblem}, Vern9{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, OrdinaryDiffEqCore.InterpolationData{ODEFunction{true, SciMLBase.AutoSpecialize, typeof(diffEquations!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, Vector{Vector{PrecisionCarrier{Float64}}}, Vector{PrecisionCarrier{Float64}}, Vector{Vector{Vector{PrecisionCarrier{Float64}}}}, Nothing, OrdinaryDiffEqVerner.Vern9Cache{Vector{PrecisionCarrier{Float64}}, Vector{PrecisionCarrier{Float64}}, Vector{Float64}, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Nothing}, SciMLBase.DEStats, Nothing, Nothing, Nothing, Nothing}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(diffEquations!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, OrdinaryDiffEqVerner.Vern9Cache{Vector{PrecisionCarrier{Float64}}, Vector{PrecisionCarrier{Float64}}, Vector{Float64}, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, OrdinaryDiffEqCore.DEOptions{PrecisionCarrier{Float64}, PrecisionCarrier{Float64}, PrecisionCarrier{Float64}, PrecisionCarrier{Float64}, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{PrecisionCarrier{Float64}, DataStructures.FasterForward}, DataStructures.BinaryHeap{PrecisionCarrier{Float64}, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, LinRange{PrecisionCarrier{Float64}, Int64}, Tuple{}}, Nothing, Float64, Nothing, OrdinaryDiffEqCore.DefaultInit, Nothing})
    @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/zs1s7/src/solve.jl:620
  [8] #__solve#62
    @ ~/.julia/packages/OrdinaryDiffEqCore/zs1s7/src/solve.jl:7 [inlined]
  [9] __solve
    @ ~/.julia/packages/OrdinaryDiffEqCore/zs1s7/src/solve.jl:1 [inlined]
 [10] solve_call(_prob::ODEProblem{Vector{PrecisionCarrier{Float64}}, Tuple{PrecisionCarrier{Float64}, PrecisionCarrier{Float64}}, true, Tuple{Vector{PrecisionCarrier{Float64}}, Int64, Vector{PrecisionCarrier{Float64}}, Vector{PrecisionCarrier{Float64}}}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(diffEquations!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, @Kwargs{saveat::LinRange{PrecisionCarrier{Float64}, Int64}}, SciMLBase.StandardODEProblem}, args::Vern9{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{abstol::PrecisionCarrier{Float64}, reltol::PrecisionCarrier{Float64}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/yhgdI/src/solve.jl:667
 [11] solve_call
    @ ~/.julia/packages/DiffEqBase/yhgdI/src/solve.jl:624 [inlined]
 [12] #solve_up#45
    @ ~/.julia/packages/DiffEqBase/yhgdI/src/solve.jl:1199 [inlined]
 [13] solve_up
    @ ~/.julia/packages/DiffEqBase/yhgdI/src/solve.jl:1177 [inlined]
 [14] #solve#43
    @ ~/.julia/packages/DiffEqBase/yhgdI/src/solve.jl:1089 [inlined]
 [15] solve_equations()
    @ Main ~/repos/PrecisionCarriers.jl/notebooks/jl_notebook_cell_df34fa98e69747e1a8f8a730347b8e2f_W0sZmlsZQ==.jl:16
 [16] top-level scope
    @ ~/repos/PrecisionCarriers.jl/notebooks/jl_notebook_cell_df34fa98e69747e1a8f8a730347b8e2f_W0sZmlsZQ==.jl:64

```

I realize that `convert` should always return the requested type, so I patched that behavior out of my package now. With the “normal” `convert` behavior, it now runs but reports basically zero epsilons everywhere, which I find hard to believe. Also, it still apparently saves a `Vector{Float64}` somewhere.

You can get the crashing behavior again by changing the first convert overload in `src/functionality/conversion.jl` to

```julia
function Base.convert(::Type{T}, p::P) where {T <: AbstractFloat}
    return P{T}(convert(T, p.x), p.big)
end

```