### Author:  [@yoshi](https://discourse.julialang.org/u/yoshi)
#### Post date: [August 2, 2023, 6:48pm UTC](https://discourse.julialang.org/t/callback-on-internal-stages-of-tsit5/102415/1 "2023-08-02T18:48:53Z")
When I step through my simulation, I can see that there are some internal stages in TSIT5 are giving values that are too large. This is to say – the variable is constrained to give values between 0 and 1 but the internal stage is giving something greater than 3. I want to create a callback function that forces recomputation of the internal stages with half the stepsize if this happens. How do I implement such a callback function?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [August 2, 2023, 9:07pm UTC](https://discourse.julialang.org/t/callback-on-internal-stages-of-tsit5/102415/2 "2023-08-02T21:07:18Z")
[https://docs.sciml.ai/OrdinaryDiffEq/stable/nonstiff/explicitrk/#OrdinaryDiffEq.Tsit5](https://docs.sciml.ai/OrdinaryDiffEq/stable/nonstiff/explicitrk/#OrdinaryDiffEq.Tsit5)
This is what the stage limiters are for.
---
### Author:  [@yoshi](https://discourse.julialang.org/u/yoshi)
#### Post date: [August 2, 2023, 9:35pm UTC](https://discourse.julialang.org/t/callback-on-internal-stages-of-tsit5/102415/3 "2023-08-02T21:35:20Z")
oye, this is a bit terse – I’m a newbie.
Is there an example somewhere of how to implement a stage\_limiter? My guess is:
```julia
using DifferentialEquations
# Define your differential equation and initial condition
function f(dx, x, p, t)
dx[1] = ...
end
x0 = [initial_value]
tspan = (t_initial, t_final)
prob = ODEProblem(f, x0, tspan)
# Define a stage_limiter function
function my_stage_limiter(u, t, integrator)
TestVals = [el[1,2] for el in integrator.k];
if sum(TestVals .< 1) != 7
integrator.dtpropose = integrator.dtpropose/2
end
return integrator
end
solver = Tsit5(stage_limiter=my_stage_limiter)
sol = solve(prob, solver)
```
I think this is wrong though because this doesn’t force recomputing the current step, it just halves the value for the next step.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [August 2, 2023, 9:38pm UTC](https://discourse.julialang.org/t/callback-on-internal-stages-of-tsit5/102415/4 "2023-08-02T21:38:29Z")
Stage limiters allow for snapping the value into the domain while keeping order. So you’d for example `integrator.u[1] > 0 && integrator.u[1] = 0` and stuff like that.
> [@yoshi](#):
>
> oye, this is a bit terse – I’m a newbie.
It’s not you it’s me. And by me I mean that we have a PR in progress to overhaul these docs.
>
Complete docs of all OrdinaryDiffEq bells and whistles should come rather soon. It’s a process to get it all in there because we expose almost everything we can, and some of that is obscure 😅.
---
### Author:  [@yoshi](https://discourse.julialang.org/u/yoshi)
#### Post date: [August 2, 2023, 9:47pm UTC](https://discourse.julialang.org/t/callback-on-internal-stages-of-tsit5/102415/5 "2023-08-02T21:47:28Z")
haha, np np! 😃 🫠
okay so, suppose I have 3 variables (variables `u[2]`, `u[3]`, and `u[4]`) I know are between 0 and 1 – then the function should read:
```julia
function my_stage_limiter(u, t, integ)
return integ.u[2] > 0 && integ.u[3] > 0 && integ.u[4] > 0 && integ.u[2] < 1 && integ.u[3] > 1 && integ.u[4] < 1
end
```
If this evaluates false during any of the internal stages, the step will recompute with a smaller stepsize?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [August 2, 2023, 9:53pm UTC](https://discourse.julialang.org/t/callback-on-internal-stages-of-tsit5/102415/6 "2023-08-02T21:53:29Z")
> [@yoshi](#):
>
> If this evaluates false during any of the internal stages, the step will recompute with a smaller stepsize?
No, stage limiters are for modifying the stage value. We don’t have a mechanism for forcing no violations between stages. @ranocha you had an algorithm you pushed for that though? Maybe we should look into adding something for this case.
---
### Author:  [@yoshi](https://discourse.julialang.org/u/yoshi)
#### Post date: [August 2, 2023, 11:46pm UTC](https://discourse.julialang.org/t/callback-on-internal-stages-of-tsit5/102415/7 "2023-08-02T23:46:22Z")
Is there a way to let all the stage values compute, then check them all for violations. If a violation exists, recompute the stage values for half the stepsize?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [August 2, 2023, 11:58pm UTC](https://discourse.julialang.org/t/callback-on-internal-stages-of-tsit5/102415/8 "2023-08-02T23:58:53Z")
The keyword argument `isoutofdomain` can take a function to reject steps by what comes out after the last stage.
---
### Author:  [@yoshi](https://discourse.julialang.org/u/yoshi)
#### Post date: [August 3, 2023, 12:08am UTC](https://discourse.julialang.org/t/callback-on-internal-stages-of-tsit5/102415/9 "2023-08-03T00:08:48Z")
Okay, so it takes a function that returns true or false. If it returns false, the step is recomputed with a smaller stepsize?
---
### Author:  [@yoshi](https://discourse.julialang.org/u/yoshi)
#### Post date: [August 3, 2023, 1:37am UTC](https://discourse.julialang.org/t/callback-on-internal-stages-of-tsit5/102415/10 "2023-08-03T01:37:11Z")
nvm, its in the documentation lol
---
### Author:  [@yoshi](https://discourse.julialang.org/u/yoshi)
#### Post date: [August 3, 2023, 1:38am UTC](https://discourse.julialang.org/t/callback-on-internal-stages-of-tsit5/102415/11 "2023-08-03T01:38:57Z")
wait, but this takes only `u,p,t` as variables. Can I also pass the integrator into the function? I need to check the internal stages.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [August 3, 2023, 2:48am UTC](https://discourse.julialang.org/t/callback-on-internal-stages-of-tsit5/102415/12 "2023-08-03T02:48:50Z")
That’s a feature we can consider. It would be a breaking change but it does sound like a reasonable extension.
---
### Author:  [@ranocha](https://discourse.julialang.org/u/ranocha)
#### Post date: [August 3, 2023, 4:34am UTC](https://discourse.julialang.org/t/callback-on-internal-stages-of-tsit5/102415/13 "2023-08-03T04:34:00Z")
Yes, we have a stage callback (how we call stage limiters in Trixi.jl) to enforce such bounds for our discretizations of conservation laws, see [Callbacks · Trixi.jl](https://trixi-framework.github.io/Trixi.jl/stable/callbacks/#Stage-callbacks).
> [@ChrisRackauckas](#):
>
> That’s a feature we can consider. It would be a breaking change but it does sound like a reasonable extension.
Sounds reasonable to me, too.
---
### Author:  [@yoshi](https://discourse.julialang.org/u/yoshi)
#### Post date: [August 3, 2023, 2:52pm UTC](https://discourse.julialang.org/t/callback-on-internal-stages-of-tsit5/102415/14 "2023-08-03T14:52:20Z")
I know minimal working examples are a thing – but if you copy and paste this big block then look at the internal stage, you’ll see what I’m trying to debug. I tried breaking up the model for the last couple weeks, but I’m introducing new errors.
```julia
######################################################################################
using DifferentialEquations, NaNMath
endd = 20000000.0
#####################################################################
### Parameters
## Sensor Parameters
alphaOffset = .045;
GF = 53;
GS = 3;
GD = 1.0;
gamma = 10e-7;
delta = 10e-4;
sensorFuseDelta = .001;
## Targets
FBar = .25;
SBar = .03;
DBar = .02;
## Homeostatic Time Scales
tauG = 4000.0;
tauHalfac = 2100.0;
tauS = 2000.0;
tauAlpha = 2000.0;
## Reversal Potentials
EK(t) = -80;
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))*NaNMath.log(CaOut/CaIn);
#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(-(((EF/(sensorFuseDelta*12))^8 + (ES/(sensorFuseDelta*1.5))^8 + (ED/(sensorFuseDelta*.2))^8)^(1/8)));
############################################################################################
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,-.001)-u[40])/tauAlpha
end
############################################################################################
#Start
ICs = [-29.337327051014498 0.387305560452983 0.520851480339422 0.15123291866400054 0.7788086343207855 0.18292525224066983 0.40518167654812637 0.40305741894377733 0.09880328448541692 0.1192742729813851 0.17090849597998528 0.4853554113680752 2.455375230646913 79.93706465236964 0.2873998785854699 0.0318308211031974 0.15906997203306664 99.67308490480416 0.8899609102246536 2.6818203693916858 -6.016940524541976 6.0497029078690145 -0.08934230228658809 0.1845828676421055 0.11586968152835098 0.7364356646257856 2.4826107420682773 -6.238217364811495 -1.966199265610444 -2.1624666424135977 2.014861831075435 4.143056495983624e-7 0.999901341448654 0.16812145388254565 0.5889709108552483 0.4212141706080315 0.24284481424824048 -0.03505188961118154 -0.16940390920534412 1.0];
tspan = (0.0,endd);
prob = ODEProblem(f,ICs,tspan);
function dmnCheck(u,p,t)
return u[1] > 200.0 || u[1] < -200.0
end
integ = init(prob,Tsit5(),maxiters=1e10,isoutofdomain=dmnCheck);
solve!(integ)
```
You’ll see that `isoutofdomain` isn’t catching the failure. If you look at `integ.k` you should see an `NaN`. Some of those other values are also whacky. For instance, `1.5e74` is kinda big and `u[2]` should be between 0 and 1.
```julia
julia> integ.k
7-element Vector{Matrix{Float64}}:
[369.6253619957377 4.780405955884058 … 6.00843666970291e-6 0.0]
[558.5762304399766 3.6783055725425826 … 6.01891607017653e-6 0.0]
[163.5732051281901 2.7015515288225345 … 6.02970035701271e-6 0.0]
[19068.06445943883 -0.7991962924634176 … 6.066390790152989e-6 0.0]
[-6.30053467595449e17 -0.8521960049675845 … 6.086972935145541e-6 0.0]
[1.584810201744311e74 -2.2604713798417015 … 6.089445492965473e-6 0.0]
[NaN 13.454079037895752 … 6.420614871497879e-6 0.0]
```
Is there some workaround I can use right now that’s apparent to y’all? Or suggest what to modify directly in the code base so I can hack this thing and get it working?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [August 3, 2023, 3:08pm UTC](https://discourse.julialang.org/t/callback-on-internal-stages-of-tsit5/102415/15 "2023-08-03T15:08:04Z")
The k are the derivatives, not the stage values.
---
### Author:  [@yoshi](https://discourse.julialang.org/u/yoshi)
#### Post date: [August 3, 2023, 3:24pm UTC](https://discourse.julialang.org/t/callback-on-internal-stages-of-tsit5/102415/16 "2023-08-03T15:24:36Z")
Wait really? Is there some documentation somewhere on this stuff? I’ve just been guessing what the outputs of `fieldnames(typeof(integ))` to try to debug.
---
### Author:  [@TheLateKronos](https://discourse.julialang.org/u/TheLateKronos)
#### Post date: [August 3, 2023, 3:39pm UTC](https://discourse.julialang.org/t/callback-on-internal-stages-of-tsit5/102415/17 "2023-08-03T15:39:39Z")
> [@yoshi](#):
>
> `fieldnames(typeof(integ))`
Sidenote: I also used this for a long time time to inspect fields/properties. It turns out that the base funtion `propertynames` more or less does the same, and it feels more appropriate to me.
---
### Author:  [@yoshi](https://discourse.julialang.org/u/yoshi)
#### Post date: [August 3, 2023, 3:45pm UTC](https://discourse.julialang.org/t/callback-on-internal-stages-of-tsit5/102415/19 "2023-08-03T15:45:32Z")
Anyway, there should be bounds on derivatives too, like some of these must be nonnegative and some of these must be nonpositive. `u[2]` should not then exhibit behavior that is both positive and negative in `integ.k`.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [August 4, 2023, 12:37pm UTC](https://discourse.julialang.org/t/callback-on-internal-stages-of-tsit5/102415/20 "2023-08-04T12:37:34Z")
> [@yoshi](#):
>
> Wait really? Is there some documentation somewhere on this stuff?
There is no intention to document things that are very internal to the algorithms. You can poke around all you want since the code is open (that’s why it’s open) and interpret it to use it at your own risk (though I will say the `k` stuff hasn’t changed in about 6 years for the ERK methods so it’s pretty stable).
Note that it’s not the derivatives for all algorithms so use it at your own risk and check what it means.
Though you’re trying to accomplish something rather odd. I think a quick fix would just be to make a dispatch of a version of the method where the limiter is passed the integrator and that should be sufficient for what you need.
---
### Author:  [@yoshi](https://discourse.julialang.org/u/yoshi)
#### Post date: [August 4, 2023, 1:23pm UTC](https://discourse.julialang.org/t/callback-on-internal-stages-of-tsit5/102415/21 "2023-08-04T13:23:51Z")
I mean I tried passing a limiter to the integrator to check the u’s, it’s not working. I’m not sure what else to try.
Im trying to figure out where the NaN is computed in integ.k.
Where in the code is the tsit5 code implemented and where is integ.k loaded? Or if you don’t know off the top of your head, where would you start looking?
I looked up tsit5 as a string literal in the sciml code base, but that wasn’t helpful.
[Next page](https://discourse.julialang.org/t/callback-on-internal-stages-of-tsit5/102415.md?page=2)