Value of variable not changing.

In the following code the value mdot is set to change once time t reaches a certain value but when i run it it does not change. Could anyone explain why this is happening?

clear all;
%Projectile motion with drag and thrust solution using RK4
%variation of density with height
%x'' = FDx + Tx, y'' = -g + FDy + Ty Equations to be solved

g=9.807;
vel=1000; th_deg=60;m=80; %input 
cd=0.75; rho=1.225; S=2.176*10^(-4);
k=0.5*rho*cd*S; % constant k used in drag force F=kv^2
mdot = 0.2; %burn rate
Isp1 = 280; Isp2=480; %stage wise specific impulse
Isp=Isp1;
ms1=20; %structural mass of first stage

x0=0; y0=0; %initial condition
t0=0; tf=10000; %time span
vx=vel*cosd(th_deg); %velocity along x
vy=vel*sind(th_deg); %velocity along y

beta=deg2rad(th_deg);%angle made by drag force with horizontal

%transforming second order differential equation to first order
%x'=u represented by fx 
%u'=-(k/m)*(u^2+v^2)*cos(beta)+mdot*g*Isp*cos(beta)/m represented by fu
%y'=v represented by fy
%v'=-g-(k/m)*v*(u^2+v^2)*sin(beta) +mdot*g*Isp*sin(beta)/m represented by fv
fx=@(t,x,u) u;
fu=@(t,x,u,v)(-(k/m)*(u^2+v^2)*cos(beta)+(mdot*g*Isp*cos(beta))/m);
fy=@(t,y,v) v;
fv=@(t,y,u,v)(-g-(k/m)*(u^2+v^2)*sin(beta)+(mdot*g*Isp*sin(beta))/m);

t(1)=0;
x(1)=0;y(1)=0;
u(1)=vx;v(1)=vy;
h=0.01;
N=ceil((tf-t(1))/h);
for j=1:N
    t(j+1)=t(j)+h;
   
    k1x=fx(t(j),x(j),u(j));
    k1y=fy(t(j),y(j),v(j));
    k1u=fu(t(j),x(j),u(j),v(j));
    k1v=fv(t(j),y(j),u(j),v(j));
    
    k2x=fx(t(j)+h/2,x(j)+h/2*k1x,u(j)+h/2*k1u);
    k2y=fy(t(j)+h/2,y(j)+h/2*k1y,v(j)+h/2*k1v);
    k2u=fu(t(j)+h/2,x(j)+h/2*k1x,u(j)+h/2*k1u,v(j)+h/2*k1v);
    k2v=fv(t(j)+h/2,y(j)+h/2*k1y,u(j)+h/2*k1u,v(j)+h/2*k1v);
    
    k3x=fx(t(j)+h/2,x(j)+h/2*k2x,u(j)+h/2*k2u);
    k3y=fy(t(j)+h/2,y(j)+h/2*k2y,v(j)+h/2*k2v);
    k3u=fu(t(j)+h/2,x(j)+h/2*k2x,u(j)+h/2*k2u,v(j)+h/2*k2v);
    k3v=fv(t(j)+h/2,y(j)+h/2*k2y,u(j)+h/2*k2u,v(j)+h/2*k2v);
    
    k4x=fx(t(j)+h,x(j)+h*k3x,u(j)+h*k3u);
    k4y=fy(t(j)+h,y(j)+h*k3y,v(j)+h*k3v);
    k4u=fu(t(j)+h,x(j)+h*k3x,u(j)+h*k3u,v(j)+h*k3v);
    k4v=fv(t(j)+h,y(j)+h*k3y,u(j)+h*k3u,v(j)+h*k3v);
    
    x(j+1)=x(j)+h/6*(k1x+2*k2x+2*k3x+k4x);
    y(j+1)=y(j)+h/6*(k1y+2*k2y+2*k3y+k4y);
    u(j+1)=u(j)+h/6*(k1u+2*k2u+2*k3u+k4u);
    v(j+1)=v(j)+h/6*(k1v+2*k2v+2*k3v+k4v);
    
    m=m-mdot*h;
    
    %condition for burn time for first stage
    switch(t(j+1))
        case 120
        mdot=0;
        m=m-ms1;%first stage separation
    
        %no thrust phase
        %second stage ignition 
        case 240
        mdot=0.8;
        Isp=Isp2;
   
        %condition for burn time for second stage
        case 360
        mdot=0;
    end
    
    %Density variation with height 
    if(y(j+1)<=11000)
     T=15.04 - .00649*y(j+1);
     p=101.29*((T+273.1)/288.05)^5.256;
    end
    
    if(y(j+1)>11000 && y(j+1)<=25000) 
     T=-56.46;
     p=22.65*exp(1.73-0.000157*y(j+1));
    end
    
    if(y(j+1)>25000)
     T=-131.21+0.00299*y(j+1);
     p=2.488*((T+273.1)/216.6)^(-11.388);
    end
        
    rho=p/(0.2869*(T+273.1));
    k=0.5*rho*cd*S;
    beta = atan(v(j+1)/u(j+1));
    
    
    %condition to stop simulation when particle touches earth
    if(y(j+1)<0)
        break;
    end
end

hmax = max(y);
rmax = max(x);
plot(x,y);
grid on;
xlabel('X');
ylabel('Y');

[Sahil Kumar]

Thankyou. It was solved.

Answer 1

[Walter Roberson]

Your code assumes that if you add up enough 0.01 that the result will be an exact integer such as 120 or 240. However,

format long
t=0; for K = 1 : 12000; t = t + 0.01; end
t

t = 
     1.200000000000245e+02

t - 120

ans = 
     2.448530267429305e-11

You might be thinking this is a bug, but it is normal.

MATLAB uses IEEE 754 double precision to store (most) numbers. IEEE 754 represents numbers as a 54 bit number, divided by 2 to a power. In order to represent 1/10th exactly, there would have to be some integers (P,Q) such that P/(2^Q) = 1/10 exactly. That would require that 10*P = 2^Q for integers P, Q -- there would have to be an integer multiple of 10 that exactly equalled 2 to a power. But other than 2^0 = 1, powers of 2 must end in a non-zero digit: 2, 4, 8, 6, 2, 4, 8, 6 and so no such P, Q exist.

Which is another way of saying that double precision cannot exactly represent 1/10, for the same reason that finite decimal cannot exactly represent 1/3 . 0.3333 + 0.3333 + 0.3333 = 0.9999 and 0.3334 + 0.3334 + 0.3334 = 1.0002 so you can see that if you try to represent 1/3 as 0. 3 repeated for any finite number of decimal places, or 0 . 3 repeated then final 4, then when you add three of them you cannot get exactly 1.0 back in decimal.

[Sahil Kumar]

Thankyou. It was solved.