How to save dydt(3) equation computed inside ODE?

Hi, i have to use an ode45 with a function with some if and elseif conditions inside it. The idea is that the ODE system [dydt(1) ; dydt(2) ; dydt(3)] does change depending on the solutions y(1), y(2) and y(3) that the ode45 gives each time istant of integration. The code works the way i writed it below, but now i have to plot the dydt(3) array and i dont know how. I thought to use plot( t(1:end-1), diff(y(:,3)) ), and it also work but it gives a very bad solution (too different from the real one). I also have seen the solutions proposed here:

https://it.mathworks.com/matlabcentral/answers/464400-how-to-save-a-variable-computed-inside-ode

but it doesnt work for me and i dont know why. There is a way to save the dydt(3) that result from my function? Thanks in advance.

MATLAB Version: 9.11.0.1837725 (R2021b) Update 2 .

[t,y] = ode45(@(t,y) odefcn(t,y), [0 T_f], zeros(1,3)) ;


function dydt = odefcn(t,y)

global omega_0 N_Orbits T T_f J_0 m L sigma zeta K_d K k_theta tau_theta ...
       theta_REF theta_REF_dot T_M SAT 

  dydt = zeros(3,1) ;
 
  if y(3) == 0 
      dydt(1) = y(2);
      dydt(2) = cos(t) - (k_theta*( y(1) - theta_REF ) + k_theta*tau_theta*( y(2) - theta_REF_dot )) ;
      dydt(3) = k_theta*( y(1) - theta_REF ) + k_theta*tau_theta*( y(2) - theta_REF_dot ) ;

  elseif abs(y(3)) > SAT
   if  sign( y(3) ) == sign(   k_theta*( y(1) - theta_REF ) + k_theta*tau_theta*( y(2) - theta_REF_dot )  ) 
      dydt(1) = y(2);
      dydt(2) = exp(-t) ;
      dydt(3) = 0 ;
   elseif sign( y(3) ) ~= sign(  k_theta*( y(1) - theta_REF ) + k_theta*tau_theta*( y(2) - theta_REF_dot )  ) 
      dydt(1) = y(2);
      dydt(2) = cos(t) - (k_theta*( y(1) - theta_REF ) + k_theta*tau_theta*( y(2) - theta_REF_dot )) ;
      dydt(3) = k_theta*( y(1) - theta_REF ) + k_theta*tau_theta*( y(2) - theta_REF_dot ) ;
   end

  elseif   abs(y(3)) <= SAT
      dydt(1) = y(2);
      dydt(2) = cos(t) - (k_theta*( y(1) - theta_REF ) + k_theta*tau_theta*( y(2) - theta_REF_dot )) ;
      dydt(3) = k_theta*( y(1) - theta_REF ) + k_theta*tau_theta*( y(2) - theta_REF_dot ) ;
  end

end

[Jan]

Please mention, what "it doesnt work for me" means. I've explained two working methods in the other thread. So I cannot explain anything else here. Post your code you have tried so far and explain, which problem occurs - e.g. post a copy of the complete error messages.

By the way, ODE45 is designed to integrate smooth functions. IF commands and SIGN are strong hints for functions, which are not smooth. Then the stepsite controller of ODE45 can drive mad and the final value is not a result, but dominated by accumulated rouding errors. If such a calculation occurs in a publication of PhD thesis, it would be a severe scientific mistake. Don't do this. ODE45 discontinuities

Whenever the function to be integrated has a knick or a jump, stop the integration by an event function and restart it with the new parameters or functions.

Another hint: If you have checked "if abs(y(3)) > SAT" already, checking "elseif abs(y(3)) <= SAT" is a waste of time only. elseif is sufficient.

Global variables are a shot in your knee. Use anonymous functions instead to provide parameters for the intergration: Answers: Anonymous for params

[Alessandro Castello]

Sorry i know it is not a perfect code, i just need it to work in the way i need. I tryed to run the code below and it doesnt work:

[Walter Roberson]

Yes there is, but it is almost always the wrong thing to do. ode45 evaluates the function at a variety of locations in order to predict how the function evolves, and then it tests the prediction against a different location. If the test and predictions are too different then the step is rejected and a smaller step is used; upon success a larger step will be used. It follows that except for systems that evolve as a polynomial of low degree, that ode45 will increase the step size until there is a prediction failure, so failed steps are normal. Therefore the saved history is pretty much guaranteed to include rejected locations. This is not what you want to plot.

[Jan]

@Alessandro Castello: Yes, this was a typo in my former answer. It is fixed now:

[t, y] = ode45(@fcn, [0, 2*pi], [0, 0]);
[~, p] = fcn(t.', y.');

function [dy, p] = fcn(t, y)
p  = sin(y(2, :));   % <- The wanted parameter
dy = [2 .* t; p ./ y(1, :)];
end

This method does not have the problem Walter mentions. It is evaluated for the acepted steps only. Using the OutputFcn is valid also.

Answer 1

[Torsten]

[t,y] = ode45(@(t,y) odefcn(t,y), [0 T_f], zeros(1,3)) ;

dydt = zeros(size(y));
for i = 1:numel(t)
  dydt(i,:) = odefcn(t(i),y(i,:));
end

plot(t,dydt(:,3))