How to plot the range of a model rocket?

Right now i have this code that plots the altitude with respect to time, but i want to plot the altitude with respect to range of the rocket (i.e. distance in the x-direction vs. distance in the y-direction) by chancing the launch angle. Any help is greatly appreciated! Here's what Ive got:

clear,clc;
launch_angle = 0;
dt = 0.001;
t_delay = 6; %delay time [s]
d_v = 0.0532; % vehicle diameter [m]
A_v = (pi/4)*(d_v).^2; % vehicle area [m^2]
d_p = 0.771144; % parachute diameter [m]
A_p = (pi/4)*(d_p)^2; % parachute area [m^2]
CD_V = 0.72; % vehicle drag coefficient (Will confirm after CFD)
CD_P = 1.5; % parachute drag coefficient
g = 9.81; %[m/s^2]
density = 1.225; %[kg/m^3]
I_tot = 62.2; % total impulse [Ns]
m_prop = 0.0431; %propellant mass [kg]
w_prop = m_prop*g; % propellant weight [N]
Isp_avg = I_tot/w_prop; % average specific impulse [s]
for j = 1
    time = 0;
    mass = 0.456; % initial mass [kg]
    if thrust(time) > 0
        t_burn = 2.3; % burn time [s] *provided by manufacturer*
    end
    V_0 = 0; %velocity
    Y_0 = 0; %altitude
    i = 1;
    
    cont = true;
    while cont
        if time == 0
            Y(i) = Y_0;
            V(i) = V_0;
        end
        % burn phase
        if thrust(time)>0
            mdot = thrust(time)/(9.81*Isp_avg); %propellant mass flow rate [kg/s]
            mass = 0.456 - mdot*dt; %updating mass
            D = CD_V(j)*0.5*density*(V(i)^2)*A_v; % drag [N]
            dv = (thrust(time)/mass)*dt - (D/mass)*dt - 9.81.*cosd(launch_angle).*dt;
            V(i + 1) = V(i) + dv; 
            Y(i + 1) = Y(i) + V(i)*dt;
        end
        % coast phase
        if thrust(time)<=0 && (time < t_delay+t_burn)
            D = CD_V(j)*0.5*density*(V(i)^2)*A_v;
            if(V(i) > 0 )
                D = -D;
            end
            mass = 0.456-m_prop; % initial mass - prop mass [kg]
            dv = (D/mass)*dt - 9.81.*cosd(launch_angle).*dt;
            V(i + 1) = V(i) + dv; 
            Y(i + 1) = Y(i) + V(i)*dt;
        end
        % descent phase
        if (time > t_delay+t_burn)
            D = CD_P(j)*0.5*density*(V(i)^2)*A_p;
            if(V(i) > 0 )
                D = -D;
            end
            dv = (D/mass)*dt - 9.81.*cosd(launch_angle).*dt;
            V(i + 1) = V(i) + dv; 
            Y(i + 1) = Y(i) + V(i + 1)*dt;
        end
        if(Y(i)<0)
            cont = false;
        end
        i = i + 1;
        time = time + dt;
    end
    t = 0:dt:(i - 1)*dt;
end
figure(1)
sgtitle ('Simulated Trajectory')
set(gcf,'position',[500 200 1000 600])
plot(t,Y)
set(gca,'XMinorTick','on','YMinorTick','on')
grid on,xlabel('Time [s]'),xlim([0 140]),ylim([0 500]),...
    ylabel('Altitude [m]'),xline(t_burn,':','burn time'),...
    xline(t_delay+t_burn,':','delay time'),xline(time,':','landing time')
% Data on AeroTech F26-6FJ
function T = thrust(time)
temp = [ 0.041	38.289
    0.114	36.318
    0.293	34.347
    0.497	32.939
    0.774	32.376
    1	31.25
    1.254	28.716
    1.498	25.338
    1.743	22.241
    2.003	17.737
    2.077	15.484
    2.304	5.349
    2.484	1.689
    2.61	0];
T = interp1(temp(:,1),temp(:,2),time,'linear','extrap'); 
end 

Answer 1

[Parsa]

Hi Garrett

According to your code, I didnt find any expression for velocity along x-direction. I think, you should define such an expression, with respect to initial horizontal position and Vx(t) for any integration time step.

I also recommend take a look at the link below, to find a comprehensive analysis of rocket motion, which contains a great Matlab code for that. Hope it will be useful.

MATLAB Rocket Trajectory Simulation (kevinkparsons.com)