clear
h0=39*1000;
tmin=0;

figure(1)
[t,y]=ode45('gravf',[tmin 90],[h0;0]);  % Const g
plot(t,y(:,1)/1000,'k')
hold on

[t,y]=ode45('dgravf2',[tmin 500],[h0;0]);% General rho
plot(t,y(:,1)/1000,'k:')
t2=linspace(0,300,25);
for i=1:25
    y2(i)=1.585;
end
plot(t2,y2,'k-.')
axis([0,300,0,40])
xlabel('t (s)')
ylabel('h (km)')
text(238, 3,['t=238.8 s']);
legend('No drag','General Drag','End Freefall','location','northeast');
hold off

figure(2)
plot(t,y(:,2),'k',238.8,-64.01,'ro')
xlabel('t (s)')
ylabel('v (m/s)')
title('Velocity vs Time')
axis([0,300,-400,0])


figure(3)
plot(y(:,1)/1000,y(:,2),'k',1.585,-64.01,'ro')
xlabel('h (km)')
ylabel('v (m/s)')
title('Velocity vs Altitude')
axis([0,40,-400,0])