How to fix timestep in this code?

Hello everyone, I am new to Matlab and working on a code...

I have written this code.

dt = 1; %time step
tmax = 200; %max time
t = 0:dt:tmax;
G = (6.67408*10^-11); %Gravitaionals constant m^3.kg^-1.s^-2
Mearth = (5.9722*10^24); %Mass of earth kg
R = 6371; %Radius of the earth 6371 km
A = 75; %Area of the rocket m^2
Cd = 0.4; %Drag coefficient for rocket 
h = 0; %Initial altitude
Mempty = 54000; %Mass of the rocket Empty kg
Minitial = 894000; %Mass of the rocket Full (With propellant) kg
Ve = 4500; %Exhaust gas speed  ms^-1
dm = 5000; %Rate of change of mass kg/s (constant)
V = zeros(1,length(t));
V(1) = 1;

 % using dx/dt ~= (x(i)-x(i-1))/dt Discretization example
 % speed(i) = speed(i-1) + time_step*force/mass; Remember
 
 m = zeros(1,length(t));
 m(1) = Minitial; %Mass of the rocket kg
 h = zeros(1,length(t));
 h(1) = 1 %Altitude
 p = zeros(1,length(t));
 p(1) = 1.225; %Density of air sea lvl kg/m^3

for n=2:length(t)
    if m(n-1) >= Mempty;
        m(n) = Minitial-(dt*dm*n);
    else m(n-1) = Mempty;
    end
    V(n) = V(n-1)+dt*(Ve*dm-G*((Mearth*m(n-1))/(h(n-1))^2)-0.5*p(n-1)*A*(V(n-1))^2*Cd)/m(n-1);
    h(n) = V(n)*n;
    p(n) = 1.225*10^(-3*h(n-1)/50000);
end     
plot(t,V)

    %V(n) = V(n-1)+dt*(Ve*dm-G*((Mearth*m(n))/r^2)-0.5*p*A*V(n-1)^2*Cd
    %V(n)=V(n-dt)*(1+t(n-1)*dt);
 
 

This code is supposed to calculate the velocity and altitude of a rocket flying straight up. However, for some reason as it got evident from the graph my V(n) values seem to be calculate for only 2 times steps and somehow ignores tmax variable? My assumption is that there is something wrong in my V(n) Equation, line 32.

For my V(n) I am using this

and the graph that I am getting is

which is clearly not what I am trying to get here.

I would really appreciate if someone could help me out to fix the V(n) in order to find the velocity.

Answer 1

[Gifari Zulkarnaen]

You can check your equation part by part to check whether your code or calculation is correct. For example, for n=2:

Ve*dm = 22500000 --> not weird (I don't know your field, so I don't know whether it is correct)

G*((Mearth*m(n-1))/(h(n-1))^2) = 3.5634e+20 --> Suspicious, is it correct? Maybe the r^2 in your formula means (radius_earth + h_rocket)^2...

And so on, check it part by part