Nested Function is not working correctly with ode45

I am trying to use a nested function which uses multiple if conditions to represent time periods and then I integrated using ode45 but it results in a incorrect graph I have no idea why. I am aiming to integrate between 4 time periods one between 0 and tb1, tb1 to tb+tc, tc+tb1 to tb2+tc+tb1 and greater than tb2+tc+t

This is my code:

function EGA215_1234567_A1_Q2
m01=600;
mf1=350;
T1=6500;
Isp1=270;
m02=100;
n2=1.6;
me2=0.5;
Isp2=290;
g0=9.81;
T2=Isp2*g0*me2;
tb1=101.87;
tc=10;
tb2=75;
%boost stage 1 
c1=Isp1*g0;
me1=T1/c1; 
tspan1=[0 300];
v0=0;
[t,v]=ode45(@rocket_eqn,tspan1,v0); 
figure(1)
plot(t,v)
    function vdot = rocket_eqn(t, v)
        if t<=tb1
            vdot=6500/(m01-me1*t)-g0;
        elseif (tb1<=t)<=tb1+tc
            vdot=-g0;
        elseif (tb1+tc<=t)<=tb1+tc+tb2
            vdot=T2/(m02-me2*t)-g0;  
        elseif tb2+tc+tb1<=t 
            vdot=-g0;
            
        end
    end 
end

Answer 1

[Walter Roberson]

In MATLAB, (A<=B)<=C does not mean to test whether B is between A and C. Instead, it means to test whether A<=B and if so emit 1 and otherwise emit 0, and then to compare that 0 or 1 to C.

EGA215_1234567_A1_Q2

function EGA215_1234567_A1_Q2
m01=600;
mf1=350;
T1=6500;
Isp1=270;
m02=100;
n2=1.6;
me2=0.5;
Isp2=290;
g0=9.81;
T2=Isp2*g0*me2;
tb1=101.87;
tc=10;
tb2=75;
%boost stage 1 
c1=Isp1*g0;
me1=T1/c1; 
tspan1=[0 300];
v0=0;
[t,v]=ode45(@rocket_eqn,tspan1,v0); 
figure(1)
plot(t,v)
    function vdot = rocket_eqn(t, v)
        if t<=tb1
            vdot=6500/(m01-me1*t)-g0;
        elseif t<=tb1+tc
            vdot=-g0;
        elseif t<=tb1+tc+tb2
            vdot=T2/(m02-me2*t)-g0;  
        elseif tb2+tc+tb1<=t 
            vdot=-g0;
        else
            vdot = nan;
        end
    end 
end

[Walter Roberson]

However, when you use ode45, the mathematics of Runge-Kutta requires that the first two derivatives of your function are continuous. In the great majority of cases, when people use if inside an ode function, they have not carefully ensured that the first and second derivatives are continuous at the boundaries.

If you are lucky, when you use if inside of an ode function, you will receive an error message indicating that ode45 was unable to integrate over a singularity.

If you are not lucky, then instead you will simply get the wrong result and not realize that it is the wrong result.

The above plot is not the right result.

You should study the ballode example to see how to use event functions.

Answer 2

[Torsten]

You can also get an analytical solution for v since your piecewise equation is easily integrated.

EGA215_1234567_A1_Q2

function EGA215_1234567_A1_Q2
m01=600;
mf1=350;
T1=6500;
Isp1=270;
m02=100;
n2=1.6;
me2=0.5;
Isp2=290;
g0=9.81;
T2=Isp2*g0*me2;
tb1=101.87;
tc=10;
tb2=75;
%boost stage 1 
c1=Isp1*g0;
me1=T1/c1; 
%tspan1=[0 300];
v0=0;
tspan1 = [0 tb1];
iflag = 1;
[t1,v1]=ode45(@(t,y)rocket_eqn(t,y,iflag),tspan1,v0); 
tspan2 = [tb1 tb1+tc];
iflag = 2;
[t2,v2]=ode45(@(t,y)rocket_eqn(t,y,iflag),tspan2,v1(end));
tspan3 = [tb1+tc tb1+tc+tb2];
iflag = 3;
[t3,v3]=ode45(@(t,y)rocket_eqn(t,y,iflag),tspan3,v2(end));
tspan4 = [tb1+tc+tb2 300];
iflag = 4;
[t4,v4]=ode45(@(t,y)rocket_eqn(t,y,iflag),tspan4,v3(end));
t = [t1;t2;t3;t4];
v = [v1;v2;v3;v4];
figure(1)
plot(t,v)
    function vdot = rocket_eqn(t, v,iflag)
        if iflag==1
            vdot=6500/(m01-me1*t)-g0;
        elseif iflag==2
            vdot=-g0;
        elseif iflag==3
            vdot=T2/(m02-me2*t)-g0;  
        elseif iflag==4 
            vdot=-g0;
        end
    end 
end