# solve_ivp method=ODE23 time step not decreasing in order

My time step with the function scipy.integrate.solve_ivp is not decreasing in t_span fluctuating (reaching values below or higher than t_end) which is a problem for my application. Can anyone tell me if, first, it is a normal situation and, second, if it is possible to make the time step continuously decrease in t_span?

Here is an example of the time step that the function is passing through

-32.0
-24.572060540292647 <--- This first drop is the problem
-31.895622059576212
-31.843433089364318
-31.791244119152424
-31.395622059576212
-31.197811029788106
-31.0
t_span = [-32, -31]


Here the part of my code using this function. I use THETA to compute X and depending on the value of X, I know if I am in the compression stroke, combustion phase, or expansion stroke but this time drop put me in the combustion phase but I am actually in the compression stroke.

def rates( THETA, Y):
YPRIME = zeros(NY);
VOL, X , EM = auxiliary(THETA);
print(THETA)
M = EM*MNOT;
VOL, DV, Ac = volume(THETA)
AA = (DV + VOL*BLOWBY/OMEGA)/M;
C1 = HEAT*(Ac)/OMEGA/M/1000;
C0 = sqrt(X);
P = Y[0];
TB = Y[1];
TU = Y[2];

if ( X > 0.999 ):
#  EXPANSION
ierr, YB, HL, xxx, xxx, VB, xxx, CP, xxx, DVDT, DVDP = ecp( TB, P, PHI, fuel_type );
if ( ierr != 0 ):
print('Error in ECP({0}, {1}, {2}): {3}\n'.format(TB, P, PHI, ierr))

BB = C1/CP*DVDT*TB*(1-TW/TB);
CC = 0;
DD = 1/CP*TB*DVDT**2 + DVDP;
EE = 0;

YPRIME[0] = (AA + BB + CC)/(DD + EE);
YPRIME[1] = -C1/CP*(TB-TW) + 1/CP*DVDT*TB*YPRIME[0];
YPRIME[2] = 0;

elif ( X > 0.001 ):
#  COMBUSTION
xxx, HU, xxx, xxx, VU, xxx, CPU, xxx, DVDTU, DVDPU = farg( TU, P, PHI, F, fuel_type );
ierr, YB, HB, xxx, xxx, VB, xxx, CPB, xxx, DVDTB, DVDPB = ecp( TB, P, PHI, fuel_type );
if ( ierr != 0 ):
print('Error in ECP({0}, {1}, {2}): {3}\n'.format(TB, P, PHI, ierr))
BB = C1*(1/CPB*TB*DVDTB*C0*(1-TW/TB) + 1/CPU*TU*DVDTU*(1-C0)*(1-TW/TU));
DX = xb(THETA,THETAS,THETAB) [1]
CC = -(VB-VU)*DX - DVDTB*(HU-HB)/CPB*(DX-(X-X**2)*BLOWBY/OMEGA);
DD = X*(VB**2/CPB/TB*(TB/VB*DVDTB)**2 + DVDPB);
EE = (1-X)*(1/CPU*TU*DVDTU**2 + DVDPU);
HL = (1-X**2)*HU + X**2*HB;

YPRIME[0] = (AA + BB + CC)/(DD + EE);
YPRIME[1] = -C1/CPB/C0*(TB-TW) + 1/CPB*TB*DVDTB*YPRIME[0] + (HU-HB)/CPB*(DX/X - (1-X)*BLOWBY/OMEGA);
YPRIME[2] = -C1/CPU/(1+C0)*(TU-TW) + 1/CPU*TU*DVDTU*YPRIME[0];

else:
#  COMPRESSION
xxx, HL, xxx, xxx, xxx, xxx, CP, xxx, DVDT, DVDP = farg( TU, P, PHI, F, fuel_type )

BB = C1*1/CP*TU*DVDT*(1-TW/TU)
CC = 0
DD = 0
EE = 1/CP*TU*DVDT**2 + DVDP
YPRIME[0] = ( AA + BB + CC )/(DD + EE)
YPRIME[1] = 0
YPRIME[2] = -C1/CP*(TU-TW) + 1/CP*TU*DVDT*YPRIME[0]

# 1/omega is s/rad, so convert to s/deg
for JJ in range(NY) :
YPRIME[JJ] = YPRIME[JJ]*pi/180

return YPRIME

def integrate( THETA, THETAE, Y, X):
sol = scipy.integrate.solve_ivp(rates, array([THETA, THETAE]), Y, method='RK23', args=(X,))
return sol


Scipy's solve_ivp has indeed an issue where it may sometimes go beyond the final time at the first step. Since you do not provide an initial step size (via the first_step argument), it attempts computing one based on a simple method (described e.g. in Hairer ang Wanner, "Solving Ordinary Differential Equations I", sect II.4). See the source code here: https://github.com/scipy/scipy/blob/f43f724827c9e2d40a04b7ab4bd2d46a2632c241/scipy/integrate/_ivp/common.py#L68