# Troubleshooting Integration

**URL:** <https://discourse.julialang.org/t/troubleshooting-integration/98957>\
**Category:** New to Julia\
**Created:** [May 16, 2023, 9:25pm UTC](https://discourse.julialang.org/t/troubleshooting-integration/98957 "2023-05-16T21:25:46Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![yoshi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yoshi/32/38259_2.png) [@yoshi](https://discourse.julialang.org/u/yoshi)\
**Post date:** [May 16, 2023, 9:25pm UTC](https://discourse.julialang.org/t/troubleshooting-integration/98957/1 "2023-05-16T21:25:46Z")

</div>

Hi,

I’m having trouble thinking about how to troubleshoot the following equation. I’ve isolated the time step that produces the error and have examined many values that the error suggests to look at (for instance, the values of ECa, Ca, ICaT, ICaS, dCa/dt - equation 13, etc). Nothing looks a miss. The next step for me was to examine the values the integrator is producing. But I’m having trouble using the integrator interface. Below is the code and things I’ve tried. Do people see what’s going on? Thanks!

Code

```julia
using DifferentialEquations

endd = 99550.9 # last timestep that works

#endd = 99551.0 # first timestep that fails

function EK(t)
    if t < 0
        return -80
    elseif t >= 0 && t < 50000.0
        return -80
    elseif t >= 50000.0 && t < 150000.0
        return -80
    else
        return -80
    end
end

#####################################################################
### Parameters

## Sensor Parameters
alphaOffset = 0.075 # alpha in between
GF = 53;
GS = 3;
GD = 1.0;
gamma = 10e-7;
delta = 10e-4;
sensorFuseDelta = .005;

## Targets
FBar = .25;
SBar = .03; 
DBar = .02;

## Homeostatic Time Scales
tauG = 4000.0;
tauHalfac = Inf;
tauS = 2000.0;
tauAlpha = 2000.0;

## Reversal Potentials
EH = -20.0;
ENa = 30.0;
EL = -50.0;

#####################################################################

### Define generic functions
## Generic Activation/Inactication Curves
sigmd(x,translate,center,reactance) = 1/(1+exp(((x-translate) + center)/reactance));
## Define the time constant function for ion channels
tauX(Volt, CT, DT, AT, BT, halfAc) = CT - DT/(1 + exp(((Volt-halfAc) + AT)/BT));
## Some intrinsic currents require a special time constant function
spectau(Volt, CT, DT, AT, BT, AT2, BT2, halfAc) = CT + DT/(exp(((Volt-halfAc) + AT)/BT) + exp(((Volt-halfAc) + AT2)/BT2));
## Define the ionic currents; q is the exponent of the activation variable m
iIonic(g, m, h, q, Volt, Erev) = g*(m^q)*h*(Volt - Erev);

#####################################################################

# Ion Channel Acitvation/Inactivation Curves
NaMinf(V,trns) = sigmd(V,trns, 25.5, -5.29); # m^3
NaHinf(V,trns) = sigmd(V,trns, 48.9, 5.18); # h
CaTMinf(V,trns) = sigmd(V,trns, 27.1, -7.20); # m^3
CaTHinf(V,trns) = sigmd(V,trns, 32.1, 5.50); # h
CaSMinf(V,trns) = sigmd(V,trns, 33.0, -8.1); # m^3
CaSHinf(V,trns) = sigmd(V,trns, 60.0, 6.20); # h
HMinf(V,trns) = sigmd(V,trns, 70.0, 6.0); # m
KdMinf(V,trns) = sigmd(V,trns, 12.3, -11.8); # m^4
KCaMinf(V,trns,IntCa) = (IntCa/(IntCa + 3.0))*sigmd(V,trns, 28.3, -12.6); # m^4
AMinf(V,trns) = sigmd(V,trns, 27.2, -8.70); # m^3
AHinf(V,trns) = sigmd(V,trns, 56.9, 4.90); # h

# Time Constants (ms)
tauNaM(V,trns) = tauX(V, 1.32, 1.26, 120.0, -25.0, trns);
tauNaH(V,trns) = tauX(V, 0.00, -0.67, 62.9, -10.0, trns)*tauX(V, 1.50, -1.00, 34.9, 3.60, trns);
tauCaTM(V,trns) = tauX(V, 21.7, 21.3, 68.1, -20.5, trns);
tauCaTH(V,trns) = tauX(V, 105.0, 89.8, 55.0, -16.9, trns);
tauCaSM(V,trns) = spectau(V, 1.40, 7.00, 27.0, 10.0, 70.0, -13.0, trns);
tauCaSH(V,trns) = spectau(V, 60.0, 150.0, 55.0, 9.00, 65.0, -16.0, trns);
tauHM(V,trns) = tauX(V, 272.0, -1499.0, 42.2, -8.73, trns);
tauKdM(V,trns) = tauX(V, 7.20, 6.40, 28.3, -19.2, trns);
tauKCaM(V,trns) = tauX(V, 90.3, 75.1, 46.0, -22.7, trns);
tauAM(V,trns) = tauX(V, 11.6, 10.4, 32.9, -15.2, trns);
tauAH(V,trns) = tauX(V, 38.6, 29.2, 38.9, -26.5, trns);

# Calcium Reverse Potential
R = 8.314*1000.0; # Ideal Gas Constant (*10^3 to put into mV)
temp = 10.0; # Temperature in Celcius; Temperature of the Sea lawl
z = 2.0; # Valence of Calcium Ions
Far = 96485.33; # Faraday's Constant
CaOut = 3000.0; # Outer Ca Concentration (uM)
ECa(CaIn) = ((R*(273.15 + temp))/(z*Far))*log(CaOut/CaIn);

# Ionic Currents (mV / ms)
iNa(g,m,h,V) = iIonic(g, m, h, 3.0, V, ENa);
iCaT(g,m,h,V,CaIn) = iIonic(g, m, h, 3.0, V, ECa(CaIn));
iCaS(g,m,h,V,CaIn) = iIonic(g, m, h, 3.0, V, ECa(CaIn));
iH(g,m,V) = iIonic(g, m, 1.0, 1.0, V, EH);
iKd(g,m,V,t) = iIonic(g, m, 1.0, 4.0, V, EK(t));
iKCa(g,m,V,t) = iIonic(g, m, 1.0, 4.0, V, EK(t));
iA(g,m,h,V,t) = iIonic(g, m, h, 3.0, V, EK(t));
iL(V) = iIonic(.01, 1.0, 1.0, 1.0, V, EL);
iApp(t) = 0;

# Sensor Equations
F(FM,FH) = GF*FM^2*FH;
S(SM,SH) = GS*SM^2*SH;
D(DM) = GD*DM^2; 
sensorFuse(EF,ES,ED) = exp(-abs(EF)/(sensorFuseDelta*10))*exp(-abs(ES)/sensorFuseDelta)*exp(-abs(ED)/sensorFuseDelta);

############################################################################################

function f(du,u,p,t)
    du[1] = -iL(u[1])-iNa(u[14],u[2],u[3],u[1])-iCaT(u[15],u[4],u[5],u[1],u[13])-iCaS(u[16],u[6],u[7],u[1],u[13])-iH(u[17],u[8],u[1])-iKd(u[18],u[9],u[1],t)-iKCa(u[19],u[10],u[1],t)-iA(u[20],u[11],u[12],u[1],t)
    du[2] = (NaMinf(u[1],u[21]) - u[2])/tauNaM(u[1],u[21])
    du[3] = (NaHinf(u[1],u[22]) - u[3])/tauNaH(u[1],u[22])
    du[4] = (CaTMinf(u[1],u[23]) - u[4])/tauCaTM(u[1],u[23])
    du[5] = (CaTHinf(u[1],u[24]) - u[5])/tauCaTH(u[1],u[24])
    du[6] = (CaSMinf(u[1],u[25]) - u[6])/tauCaSM(u[1],u[25])
    du[7] = (CaSHinf(u[1],u[26]) - u[7])/tauCaSH(u[1],u[26])
    du[8] = (HMinf(u[1],u[27]) - u[8])/tauHM(u[1],u[27])
    du[9] = (KdMinf(u[1],u[28]) - u[9])/tauKdM(u[1],u[28])
    du[10] = (KCaMinf(u[1],u[29],u[13]) - u[10])/tauKCaM(u[1],u[29])
    du[11] = (AMinf(u[1],u[30]) - u[11])/tauAM(u[1],u[30])
    du[12] = (AHinf(u[1],u[31]) - u[12])/tauAH(u[1],u[31])
    du[13] = (-.94*(iCaT(u[15],u[4],u[5],u[1],u[13])+iCaS(u[16],u[6],u[7],u[1],u[13]))-u[13]+.05)/20
    du[14] = ((1*(FBar-F(u[32],u[33])) + 0*(SBar-S(u[34],u[35])) + 0*(DBar-D(u[36])))*u[14]-gamma*u[14]^3)*u[40]/tauG
    du[15] = ((0*(FBar-F(u[32],u[33])) + 1*(SBar-S(u[34],u[35])) + 0*(DBar-D(u[36])))*u[15]-gamma*u[15]^3)*u[40]/tauG
    du[16] = ((0*(FBar-F(u[32],u[33])) + 1*(SBar-S(u[34],u[35])) + 0*(DBar-D(u[36])))*u[16]-gamma*u[16]^3)*u[40]/tauG
    du[17] = ((0*(FBar-F(u[32],u[33])) + 1*(SBar-S(u[34],u[35])) + 1*(DBar-D(u[36])))*u[17]-gamma*u[17]^3)*u[40]/tauG
    du[18] = ((1*(FBar-F(u[32],u[33])) + -1*(SBar-S(u[34],u[35])) + 0*(DBar-D(u[36])))*u[18]-gamma*u[18]^3)*u[40]/tauG
    du[19] = ((0*(FBar-F(u[32],u[33])) + -1*(SBar-S(u[34],u[35])) + -1*(DBar-D(u[36])))*u[19]-gamma*u[19]^3)*u[40]/tauG
    du[20] = ((0*(FBar-F(u[32],u[33])) + -1*(SBar-S(u[34],u[35])) + -1*(DBar-D(u[36])))*u[20]-gamma*u[20]^3)*u[40]/tauG
    du[21] = ((-1*(FBar-F(u[32],u[33])) + 0*(SBar-S(u[34],u[35])) + 0*(DBar-D(u[36])))-delta*u[21]^3)*u[40]/tauHalfac
    du[22] = ((1*(FBar-F(u[32],u[33])) + 0*(SBar-S(u[34],u[35])) + 0*(DBar-D(u[36])))-delta*u[22]^3)*u[40]/tauHalfac
    du[23] = ((0*(FBar-F(u[32],u[33])) + -1*(SBar-S(u[34],u[35])) + 0*(DBar-D(u[36])))-delta*u[23]^3)*u[40]/tauHalfac
    du[24] = ((0*(FBar-F(u[32],u[33])) + 1*(SBar-S(u[34],u[35])) + 0*(DBar-D(u[36])))-delta*u[24]^3)*u[40]/tauHalfac
    du[25] = ((0*(FBar-F(u[32],u[33])) + -1*(SBar-S(u[34],u[35])) + 0*(DBar-D(u[36])))-delta*u[25]^3)*u[40]/tauHalfac
    du[26] = ((0*(FBar-F(u[32],u[33])) + 1*(SBar-S(u[34],u[35])) + 0*(DBar-D(u[36])))-delta*u[26]^3)*u[40]/tauHalfac
    du[27] = ((0*(FBar-F(u[32],u[33])) + -1*(SBar-S(u[34],u[35])) + -1*(DBar-D(u[36])))-delta*u[27]^3)*u[40]/tauHalfac
    du[28] = ((-1*(FBar-F(u[32],u[33])) + 1*(SBar-S(u[34],u[35])) + 0*(DBar-D(u[36])))-delta*u[28]^3)*u[40]/tauHalfac
    du[29] = ((0*(FBar-F(u[32],u[33])) + 1*(SBar-S(u[34],u[35])) + 1*(DBar-D(u[36])))-delta*u[29]^3)*u[40]/tauHalfac
    du[30] = ((0*(FBar-F(u[32],u[33])) + 1*(SBar-S(u[34],u[35])) + 1*(DBar-D(u[36])))-delta*u[30]^3)*u[40]/tauHalfac
    du[31] = ((0*(FBar-F(u[32],u[33])) + -1*(SBar-S(u[34],u[35])) + -1*(DBar-D(u[36])))-delta*u[31]^3)*u[40]/tauHalfac
    du[32] = (sigmd(iCaT(u[15],u[4],u[5],u[1],u[13])+iCaS(u[16],u[6],u[7],u[1],u[13]),0,14.8,1) - u[32])/.5
    du[33] = (sigmd(-1*(iCaT(u[15],u[4],u[5],u[1],u[13])+iCaS(u[16],u[6],u[7],u[1],u[13])),0,-9.8,1) - u[33])/1.5
    du[34] = (sigmd(iCaT(u[15],u[4],u[5],u[1],u[13])+iCaS(u[16],u[6],u[7],u[1],u[13]),0,7.2,1) - u[34])/50
    du[35] = (sigmd(-1*(iCaT(u[15],u[4],u[5],u[1],u[13])+iCaS(u[16],u[6],u[7],u[1],u[13])),0,-2.8,1) - u[35])/60
    du[36] = (sigmd(iCaT(u[15],u[4],u[5],u[1],u[13])+iCaS(u[16],u[6],u[7],u[1],u[13]),0,3,1) - u[36])/500
    du[37] = ((FBar - F(u[32],u[33])) - u[37])/tauS
    du[38] = ((SBar - S(u[34],u[35])) - u[38])/tauS
    du[39] = ((DBar - D(u[36])) - u[39])/tauS
    du[40] = (sigmd(-1*sensorFuse(u[37],u[38],u[39]),0,alphaOffset,-.01)-u[40])/tauAlpha
end

############################################################################################

ICs = [-50.0 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.01 0.8320755796233367 0.046612330545369285 0.4039904957344005 0.7335731823371116 0.516374395070418 0.7517868565780351 0.9012425121331321 0.3795887585487998 -0.49230754975054447 -0.650925587096205 0.3021172776981724 -0.2266327552455598 -0.1602097765267385 0.2278426433885099 -0.11586687638108684 -0.5353572538829641 0.46427761017041913 -0.626348202330232 0.0 1.0 0.0 1.0 0.0 1.0 1.0 1.0 0.5];

tspan = (0.0,endd);
SaveTimes = (0.0:.1:endd); #in milliseconds

prob = ODEProblem(f,ICs,tspan);
@time sim = solve(prob,Tsit5(), maxiters=1e10, saveat=SaveTimes);

```

When run for t = 99551.0 I obtain the following error:

```julia
ERROR: DomainError with -867.2167835365867:
log will only return a complex result if called with a complex argument. Try log(Complex(x)).
Stacktrace:
  [1] throw_complex_domainerror(f::Symbol, x::Float64)
    @ Base.Math ./math.jl:33
  [2] _log(x::Float64, base::Val{:ℯ}, func::Symbol)
    @ Base.Math ./special/log.jl:292
  [3] log(x::Float64)
    @ Base.Math ./special/log.jl:257
  [4] ECa(CaIn::Float64)
    @ Main ./REPL[123]:1
  [5] iCaT(g::Float64, m::Float64, h::Float64, V::Float64, CaIn::Float64)
    @ Main ./REPL[125]:1
  [6] f(du::Matrix{Float64}, u::Matrix{Float64}, p::SciMLBase.NullParameters, t::Float64)
    @ Main ./REPL[137]:2
  [7] ODEFunction
    @ ~/.julia/packages/SciMLBase/BC0Kw/src/scimlfunctions.jl:1622 [inlined]
  [8] perform_step!(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Matrix{Float64}, Nothing, Float64, SciMLBase.NullParameters, Float64, Float64, Float64, Float64, Vector{Matrix{Float64}}, ODESolution{Float64, 3, Vector{Matrix{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Matrix{Float64}}, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, 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{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Float64, Tuple{}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Tuple{}}, Matrix{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, repeat_step::Bool)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/eqC5v/src/perform_step/low_order_rk_perform_step.jl:709
  [9] perform_step!
    @ ~/.julia/packages/OrdinaryDiffEq/eqC5v/src/perform_step/low_order_rk_perform_step.jl:675 [inlined]
 [10] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Matrix{Float64}, Nothing, Float64, SciMLBase.NullParameters, Float64, Float64, Float64, Float64, Vector{Matrix{Float64}}, ODESolution{Float64, 3, Vector{Matrix{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Matrix{Float64}}, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, 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{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Float64, Tuple{}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Tuple{}}, Matrix{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/eqC5v/src/solve.jl:477
 [11] #__solve#562
    @ ~/.julia/packages/OrdinaryDiffEq/eqC5v/src/solve.jl:5 [inlined]
 [12] #solve_call#28
    @ ~/.julia/packages/DiffEqBase/KouNZ/src/solve.jl:437 [inlined]
 [13] solve_up(prob::ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Matrix{Float64}, p::SciMLBase.NullParameters, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:maxiters, :saveat), Tuple{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/KouNZ/src/solve.jl:780
 [14] #solve#33
    @ ~/.julia/packages/DiffEqBase/KouNZ/src/solve.jl:760 [inlined]
 [15] top-level scope
    @ ./timing.jl:263 [inlined]
 [16] top-level scope
    @ ./REPL[142]:0

```

To use the integrator interface I add the following to the end of my code:

```julia
integrator = init(prob,Tsit5());
step!(integrator, 99550.9, stop_at_tdt = false)

```

And I receive the following error:

```julia
ERROR: MethodError: no method matching step!(::OrdinaryDiffEq.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Matrix{Float64}, Nothing, Float64, SciMLBase.NullParameters, Float64, Float64, Float64, Float64, Vector{Matrix{Float64}}, ODESolution{Float64, 3, Vector{Matrix{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Matrix{Float64}}, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, ODEFunction{true, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, 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{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, Matrix{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, ::Float64; stop_at_tdt=false)
Closest candidates are:
  step!(::SciMLBase.DEIntegrator, ::Any) at ~/.julia/packages/SciMLBase/BC0Kw/src/integrator_interface.jl:735 got unsupported keyword argument "stop_at_tdt"
  step!(::SciMLBase.DEIntegrator, ::Any, ::Any) at ~/.julia/packages/SciMLBase/BC0Kw/src/integrator_interface.jl:735 got unsupported keyword argument "stop_at_tdt"
  step!(::OrdinaryDiffEq.ODEIntegrator) at ~/.julia/packages/OrdinaryDiffEq/eqC5v/src/iterator_interface.jl:1 got unsupported keyword argument "stop_at_tdt"
  ...
Stacktrace:
 [1] top-level scope
   @ REPL[73]:1

```

---

<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:** [May 16, 2023, 9:31pm UTC](https://discourse.julialang.org/t/troubleshooting-integration/98957/2 "2023-05-16T21:31:21Z")

</div>

> [@yoshi](#):
>
> `stop_at_tdt = false`

That’s positional, not keyword.

`step!(integrator, 99550.9, false)`

---

<div class="post-metadata">

**Author:** ![contradict](https://avatars.discourse-cdn.com/v4/letter/c/ac91a4/32.png) [@contradict](https://discourse.julialang.org/u/contradict)\
**Post date:** [May 16, 2023, 10:59pm UTC](https://discourse.julialang.org/t/troubleshooting-integration/98957/3 "2023-05-16T22:59:27Z")

</div>

I think the timestepping just gets too aggressive and takes such a big step that the argument to the log becomes negative. One nice thing to trying in this case is to replace `Base.log` with `NaNMath.log`, returning `NaN` instead of raising an exception allows the integrator to notice something has gone wrong and try a smaller timestep. This lets your simulation complete when I tried it. In addition, `AutoTsit5(Rosenbrock23())` was substantially faster than `Tsit5()`

---

<div class="post-metadata">

**Author:** ![yoshi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yoshi/32/38259_2.png) [@yoshi](https://discourse.julialang.org/u/yoshi)\
**Post date:** [May 17, 2023, 12:04am UTC](https://discourse.julialang.org/t/troubleshooting-integration/98957/4 "2023-05-17T00:04:07Z")

</div>

Thanks!

1. I’m unsure how the integrator works – I want to return something at a resolution of “.1” (As I specify in `saveat`). Does the integrator use an adaptive time step then “interpolate” to this resolution?
2. What does `AutoTsit5(Rosenbrock23()) ` do?

---

<div class="post-metadata">

**Author:** ![yoshi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yoshi/32/38259_2.png) [@yoshi](https://discourse.julialang.org/u/yoshi)\
**Post date:** [May 17, 2023, 12:08am UTC](https://discourse.julialang.org/t/troubleshooting-integration/98957/5 "2023-05-17T00:08:18Z")

</div>

Thanks!

---

<div class="post-metadata">

**Author:** ![contradict](https://avatars.discourse-cdn.com/v4/letter/c/ac91a4/32.png) [@contradict](https://discourse.julialang.org/u/contradict)\
**Post date:** [May 17, 2023, 12:11am UTC](https://discourse.julialang.org/t/troubleshooting-integration/98957/6 "2023-05-17T00:11:14Z")

</div>

> Does the integrator use an adaptive time step then “interpolate” to this resolution?

Specifying [`saveat`](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/) ensures that the solution has a point at this time, the solver is free to actually step here or interpolate as needed to achieve the tolerances in the most efficient way possible.

`AutoTsit5(Rosenbrock23()) ` is a [Composite Algorithm](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#CompositeAlgorithm) that chooses among several algorithms based on the detected stiffness of the problem.