mdot = thrust(time)/(9.81*Isp_avg);
D = CD_V(j)*0.5*density*(V(i)^2)*A_v;
dv = (thrust(time)/mass)*dt - (D/mass)*dt - 9.81.*cosd(launch_angle).*dt;
Y(i + 1) = Y(i) + V(i)*dt;
if thrust(time)<=0 && (time < t_delay+t_burn)
D = CD_V(j)*0.5*density*(V(i)^2)*A_v;
dv = (D/mass)*dt - 9.81.*cosd(launch_angle).*dt;
Y(i + 1) = Y(i) + V(i)*dt;
if (time > t_delay+t_burn)
D = CD_P(j)*0.5*density*(V(i)^2)*A_p;
dv = (D/mass)*dt - 9.81.*cosd(launch_angle).*dt;
Y(i + 1) = Y(i) + V(i + 1)*dt;
sgtitle ('Simulated Trajectory')
set(gcf,'position',[500 200 1000 600])
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')
function T = thrust(time)
T = interp1(temp(:,1),temp(:,2),time,'linear','extrap');