ode45 graph issue

I am trying to make a code which uses ode45 functions to integrate 2 equations which will allow me to produce a graphs of velocity and altuitude against time.

The code is meant to simulate a 2 stage rocket which launches with a mass of 700kg , a thrust of 6500N and burns for 101.87s with the mass reducing at a rate of 2.4540kg/s before the first stage is jettisoned leaving a mass of 100kg, afterards the rocket has a coast phase of 10 seconds after the first boost phase where the thrust is zero, after the coast phase the mass of the rocket is reduced to 100kg and the second boost stage begins for 75 seconds in which the thrust is 1422.5N and the mass is reduced by 0.5kg/s until the end of the second boost stage.

After attempting to implement ode45 into my code i keep getting negative values for velocity and altitude which i assume means theres an error in my code logic somewhere. Below i have attached an image the differential equations used.

Thank you

% Define the thrust as a function of time
T = @(t) 6500*(t<=101.81) + 0*(t>101.81 & t<=111.81) + 1422.5*(t>111.81 & t<=186.81) + 1422.5*(t>186.81 & t<=261.81) + 0*(t>261.81); % N

% Define the mass as a function of time
m = @(t) 700 - 2.4540*t - 0*(t>101.81 & t<=111.81) - 600*(t>111.81 & t<=186.81) - 0.5*(t>186.81 & t<=261.81) - 37.5*(t>261.81); % kg

% Define the initial conditions
h0 = 0; % m
v0 = 0; % m/s

% Define the time span for the integration
tspan1 = [0, 101.81]; % s (first boost phase)
tspan2 = [101.81, 111.81]; % s (coast phase)
tspan3 = [111.81, 186.81]; % s (second boost phase)


% Integrate the system of equations using ode45 for the first boost phase
f1 = @(t, y) [y(2); T(t)/m(t) - 9.81];
[t1, y1] = ode45(f1, tspan1, [h0; v0]);

% Integrate the system of equations using ode45 for the coast phase
f2 = @(t, y) [y(2); -9.81];
[t2, y2] = ode45(f2, tspan2, [y1(end,1); y1(end,2)]);

% Integrate the system of equations using ode45 for the second boost phase
f3 = @(t, y) [y(2); T(t)/m(t) - 9.81];
[t3, y3] = ode45(f3, tspan3, [y2(end,1); y2(end,2)]);

% Integrate the system of equations using ode45 for

% Concatenate the time and state vectors
t = [t1; t2(2:end); t3(2:end)];
y = [y1; y2(2:end,:); y3(2:end,:)];

% Plot the altitude and velocity as functions of time
figure;
subplot(2,1,1);
plot(t, y(:,1));
xlabel('Time (s)');
ylabel('Altitude (m)');
subplot(2,1,2);
plot(t, y(:,2));
xlabel('Time (s)');
ylabel('Velocity (m/s)');

Answer 1

[Steven Lord]

Instead of defining your integrand function piecewise like this:

% Define the thrust as a function of time
T = @(t) 6500*(t<=101.81) + 0*(t>101.81 & t<=111.81) + 1422.5*(t>111.81 & t<=186.81) + 1422.5*(t>186.81 & t<=261.81) + 0*(t>261.81); % N
% Define the mass as a function of time
m = @(t) 700 - 2.4540*t - 0*(t>101.81 & t<=111.81) - 600*(t>111.81 & t<=186.81) - 0.5*(t>186.81 & t<=261.81) - 37.5*(t>261.81); % kg

I would define one function for each of your time spans and solve using the appropriate function. It would help you avoid problems like the one shown by this plot of your m function:

figure
fplot(m, [0 200])
title("m")

I assume your "- 600*(t>111.81 & t<=186.81)" term was intended to simulate the drop of the first stage, but the first stage doesn't weigh 600 pounds.