Rocket Trajectory Script R-Work

I've writen a matlab script for a model water rocket to simulate its trajectory. I've done it using lots of if loops. Is there anyone that knows how I would change it to remove the if loops and make it a more standard matlab script.

% Constants and initial conditions
bottle_pressure_psi = 120; % psi
nozzle_diameter_mm = 16.5; % mm
water_density = 997; % kg/m^3
bottle_volume_l = 7.2; % L
water_volume_l = 2; % L
gas_constant = 287.05; % J/kgK
air_temp = 288; % K
empty_mass_g = 760; % g
chute_dia = 0.9; % m
chute_Cd = 0.75; %
wind_speed = 10;
gravity = -9.81; %
drag_coefficent = 0.7; %
bottle_diametre_mm = 100; % mm
air_density_constant = 1.293; % kg/m^3
water_thrust_effi = 0.90;
gas_thrust_effi = 0.90;
time_interval = 0.0001;

% Standard Units
bottle_pressure_initial = bottle_pressure_psi * 6894.75729; % Pa
nozzle_diameter = nozzle_diameter_mm / 1000; % m
nozzle_area = (3.14159 * (nozzle_diameter * nozzle_diameter)) / 4; % m^2
bottle_volume = bottle_volume_l / 1000; % m^3
water_volume_initial = water_volume_l / 1000; % m^3
empty_mass = empty_mass_g / 1000; % kg
bottle_diametre = bottle_diametre_mm / 1000; % m
bottle_area = (3.14159 * (bottle_diametre * bottle_diametre)) / 4; % m^2
water_volume = water_volume_initial;
bottle_pressure = bottle_pressure_initial;
air_density = 0;
drag_force = 0;
velocity_initial = 0;
height_initial = 0;
x_downrange = 0;

% Initialize variables
velocity = velocity_initial;
height = height_initial;
time = 0:time_interval:30; % Adjust the time range as needed
altitude_curve = zeros(size(time));
thrust_curve = zeros(size(time));
x_curve = zeros(size(time));

% Simulation loop
for i = 1:length(time)
    if water_volume > 0 % Water Thrust Stage
        % Mass Flow Rate of Water
        mass_flow_water = nozzle_area * sqrt(2 * water_density * bottle_pressure);

        % Velocity of Water
        velocity_water = mass_flow_water / (nozzle_area * water_density);

        % Thrust of Water
        thrust = (mass_flow_water * velocity_water) * water_thrust_effi;
        thrust_total = (thrust - drag_force);

        

        % Update variables for the next iteration
        water_volume = ((water_volume * 1000) - (mass_flow_water * time_interval)) / 1000;
        water_mass = water_volume * water_density;
        total_mass = water_mass + empty_mass;
        acceleration = thrust_total / total_mass;
        velocity = velocity_initial + (acceleration * time_interval);
        velocity_initial = velocity;
        bottle_pressure = (bottle_pressure_initial * (bottle_volume - water_volume_initial)) / (bottle_volume - water_volume);
        drag_force = 0.5 * air_density_constant *drag_coefficent * bottle_area * velocity^2;

       
  
    else %Air Thrust Stage
        if air_density >= 0 % Creates and initial Air Density

            % Density of Air 
            air_density = bottle_pressure / (gas_constant * air_temp);
            mass_air_initial = air_density * bottle_volume;
        end

        % Mass Flow Rate of Air
        mass_flow_air = nozzle_area * sqrt(2 * air_density * bottle_pressure);

        % Velocity of Air
        velocity_air = mass_flow_air / (nozzle_area * air_density);

        % Thrust of Air
        thrust = (mass_flow_air * velocity_air) * gas_thrust_effi;
        thrust_total = (thrust - drag_force);

        

        % Update variables for the next iteration
        mass_air = mass_air_initial - (mass_flow_air * time_interval);
        mass_air_initial = mass_air;
        total_mass = mass_air + empty_mass;
        acceleration = (thrust_total / total_mass) + gravity;
        velocity = velocity_initial + (acceleration * time_interval);
        velocity_initial = velocity;
        air_velocity = mass_flow_air / (air_density * nozzle_area);
        air_density = mass_air / bottle_volume;
        mass_flow_air = air_density * air_velocity * nozzle_area;
        bottle_pressure = ((mass_flow_air / nozzle_area)^2) / (2 * air_density);
        drag_force = 0.5 * air_density_constant *drag_coefficent * bottle_area * velocity^2;

        


        
    end

    if velocity <= -25 && thrust  <= 1 % Parachute Altitude
        velocity = -sqrt((8 * empty_mass * (-gravity)) / ( 3.14159 * air_density_constant * chute_Cd * (chute_dia^2)));
    end

    % height calculation
    height = height + velocity * time_interval;

    x_downrange = x_downrange + wind_speed *time_interval;



    % Store altitude in the array
    altitude_curve(i) = height;

    % Store thrust in the array
    
    thrust_curve(i) = thrust;

    % Store distance downrange in the array
    x_curve(i) = x_downrange;

     % Check if thrust is zero, and break the loop if true
    if height <= 0 && thrust <= 1
        break;
    end
    
   
end



% Plotting
figure;
plot(time(1:i), thrust_curve(1:i));
title('Thrust vs Time');
xlabel('Time (s)');
ylabel('Thrust (N)');

figure;
plot(time(1:i), altitude_curve(1:i));
title('Altitude vs Time');
xlabel('Time (s)');
ylabel('Altitude (m)');

figure;
plot(time(1:i), x_curve(1:i));
title('Distance Traveled vs Time');
xlabel('Time (s)');
ylabel('Distance (m)');