How to store intermediate values of solution in ode45

I have the following ode dm/dt model. I want to extract intermediate values of t, m (the solution),u at specfic intermediate points that satisfies the condition rho_air == rho_1 and rho_air == rho_2 (It is the commented block in the following code).

My question is that how to extract these values V_star1 and V_star2 and store them so that I have access to them when I call the function in ode45? Especially how to store the solution to the ode m at that point of time? I tried adding them as an additional output to the function of myode but it did not work

Thanks

clc,clear,close all
global rho_fuel Ve Mcv Cd rho_air_0 A_r g a_max rho_1 rho_2
rho_1 = 1.0098;
rho_2 = 4.7317e-04;
a_max = 32.361945; % m/s^2
g = 9.8; % m/s^2
Cd = 0.47;
r = 3.7; % m
A_r = pi * r^2;
rho_fuel = 1000; % Kg/m^3
Ve = 3000; % m/s
rho_air_0 = 1.2754; % Kg/m^3
TurnPoint = 100000;  % m
Mbody = 60684; % Kg
Mfuel = 488370; % Kg
Mpayload = 31184; % Kg
Mcv = Mbody + Mfuel + Mpayload;
tspan = [0 60];
m0 = Mfuel;
[t,m] = ode45(@(t,m) myode(t), tspan, m0);

plot(t,m,'LineWidth',2)
legend('M_e(t)','a(t)')
xlabel('time')

function [dmdt] = myode(t,m)
global rho_fuel Ve Mcv Cd rho_air_0 A_r g a_max rho_1 rho_2
a = dudt_const(t,a_max);
a_fun = @(tt) dudt_const(tt,a_max);
u = integral( a_fun , 0 , t);

h = u * t;

H=[-1;0;1;2;3;4;5;6;7;8;9;10;11;13;15;17;20;25;30;32;35;40;45;47;50;51;60;70;71;80;84.9;89.7;100.4;105;110;120];
H=H*1000;
rho=[1.347;1.225;1.1116;1.0065;0.9091;0.8191;0.7361;0.6597;0.5895;0.5252;0.4664;0.4127;0.3639;0.2655;0.1937;0.1423;0.0880;0.0395;0.0180;0.0132;0.0082;0.0039;0.0019;0.0014;0.0010;0.00086;0.000288;0.000074;0.000064;0.000015;0.000007;0.000003;0.0000005;0.0000002;0.0000001;0.0000001];
f_rho=fit(H, rho, 'linearinterp');
rho_air=f_rho(h);


% %%%%%%%%%%%%%%%%%
% if rho_air == rho_1  
%     t_star1 = t;
%     m_star1 = m; % the solution of the ode at t = t_star1
%     u_star1 = u;
%     V_star1 = [t_star1,m_star1,u_star1];
%     V_star1
% end
% 
% if rho_air == rho_2
%     t_star2 = t;
%     m_star2 = m; % the solution of the ode at t = t_star2
%     u_star2 = u;
%     V_star2 = [t_star2,m_star2,u_star2];
%     V_star2
% end
% 
% %%%%%%%%%%%%%%%%%

Fd = 0.5 * rho_air * u^2 * A_r * Cd;
dmdt = (1/Ve) * (-Mcv * g - Fd - Mcv * a);
end
function a = dudt_const(t,a_max)    
if t == 0
    a = a_max;
end
if t > 0 
    a = a_max .* t./t;
end
end

[Jan]

Note: myode need 2 input arguments. Then:

[t,m] = ode45(@(t,m) myode(t), tspan, m0);
%                          ^ 1 argument only

Better:

[t,m] = ode45(@(t,m) myode(t, m), tspan, m0);
% Or simply:
[t,m] = ode45(@myode, tspan, m0);

Answer 1

[Walter Roberson]

Doing that would be a mistake.

Runge-Kutta methods such as ode45 always have a "current" state, and a current step size. They evaluate the ode function at a series of carefully-chosen locations "near" the current state, at distances proportionate to the current step size. The results are used to predict the value at another point. Then the ode function is evaluated at the predicted point. If the difference in value between the prediction and the calculated value is within tolerance, then the location is accepted and recorded, the step size is increased, and the location becomes the new "current" point. If the difference between prediction and actual is too great then the point is rejected and the step size is decreased and the ode* routine tries again from the previous point. This all involves evaluating at a number of points which are highly unlikely to ever be accepted themselves even if the step is successful, and certainly not if the proposal is rejected.

If you were returning values out of the ode function, you would be doing so for all of the locations that were probed, even though there was no real chance they would become current points.

Think of it sort of like walking on a curved hilly path in the dark with a narrow-beam flashlight whose batteries are running out. You shine your flashlight very briefly in careful patterns to predict which way it will be safe to walk. You know you aren't going to walk up that outcropping, but looking there and over there give you information about how steep the area is. You don't really want to keep the details of all of those locations that you glanced at: you only need the details of where you actually walked, but you needed to look at those other places too in order to figure out where to walk.

It also turns out that if you specified a tspan that has more than two entries, that the ode*() routines might happen to only ever report back about locations that it estimated based on where you did go. For example if you asked for a report every 5 seconds, then it might have picked a location at 44 seconds and another at 47 seconds, and so to report back every 5 seconds it would estimate the 45 second information based on the 44 and 47 second information.

If you only specify two elements in the time vector, then it turns out that with the default options, 3/4 of the outputs are esimates (you can tell it to only report locations that were "accepted" though.)

What you really want to know, almost all of the time, is the details associated with the locations that were reported -- locations that might never have been evaluated at, as they might be mostly estimates. So most of the time you should not try to record details while the ode*() function is running. Instead, take the locations that are output, and ask for details about those specific locations.