Why is ode45 changing the value of the state variable between calls to my dynamics function?

Hello. I am running ode45 to simulate the dynamics of a nonlinear system. Inside the dynamics function that ode45 calls, I have an if statement that sets two of the state variables (thetadot and thetaddot) to 0 if a certain condition is met. However, when ode45 calls the dynamics function again on the next time step, thetadot is NOT equal to zero. Why is the value of thetadot changing from one time step to the next?

Here is my call to ode45:

[t,y] = ode45(@(Tspan, y0) singleDynamicsCoil(Tspan,y0,mu0,coils,mag,...
    mass,w,L,h,ct,cn,Bmax), Tspan, y0, options) ; 

And here is the if statement that sets thetadot and thetaddot to 0:

if abs(Bangle - theta) < 1E-1
    disp('***ACTIVATING TH = 0 CONDITION***') ; 
    Th = 0 ; 
    Tb = 0 ; 
    thetadot = 0 ; 
    thetaddot = 0 ; 
    disp('thetadot') ; 
    disp(thetadot) ;
    disp('thetaddot') ; 
    disp(thetaddot) ; 
else
    thetaddot = (1/I) * (Tb - Th) ;
end

%%%%%%%%%%%%%%%%%%%%%%
% return accelerations 
%%%%%%%%%%%%%%%%%%%%%%
qdot = [xdot ; ydot ; thetadot ; xddot ; yddot ; thetaddot] ;   
end % this is the end of the dynamics function

I have a display statement at the top of the dynamics function here:

function qdot = singleDynamicsCoil(t,q,mu0,coils,mag,mass,w,L,h,ct,cn,Bmax) % this is the beginning of the dynamics function

    disp('****START****') ; 
    
    %%%%%%%%%%%%%%%%%
    % state variables
    %%%%%%%%%%%%%%%%% 
    x = q(1) ; 
    y = q(2) ; 
    theta = q(3)  ;
    xdot = q(4) ; 
    ydot = q(5) ; 
    thetadot = q(6) ; 
   
    disp('thetadot') ; 
    disp(thetadot) ; 

The last line, disp(thetadot) is the line that reveals that even though I set thetadot = 0 at the end of the dynamics function in the last time step, its value has been changed to some non-zero value before the next call to the dynamics function. Why is this happening?

Thank you in advance for your time and help.

[Star Strider]

How does the if block interact with the ‘singleDynamicsCoil’ function?

Also, ‘Bangle’ does not appear to be in ‘singleDynamicsCoil’, so that could also be a problem.

[Walter Roberson]

At every place of discontinuity in the equations, you must stop the ode45 call and start again from that location. See ballode example.

Answer 1

[Jim Riggs]

ode45 implements the Dormand-Prince method which is an embedded 4th / 5th order Runge-Kutta method. For every "step" in the solution, ode45 needs to execute your function six times; i.e. there are six "stages" per every solution "step". Each stage will involve a change in the time and/or the state values.

So, at the end of one step, what you have is the result of a combibnation of all six stages, not just the last call to your function.