h0=1000;
tmax=20;
tmin=0;

[t,y]=ode45('odef',[tmin .75*tmax],[h0;0]);  % Const g
plot(t,y(:,1),'k')
hold on
[t,y]=ode45('gravf',[tmin .75*tmax],[h0;0]); % Nonconst g
plot(t,y(:,1),'k.')

[t,y]=ode45('dgravf',[tmin tmax],[h0;0]);% Const rho
plot(t,y(:,1),'k--')

[t,y]=ode45('dgravf2',[tmin tmax],[h0;0]);% General rho
plot(t,y(:,1),'k:')

xlabel('t (s)')
ylabel('h (m)')
legend('Const g, C=0','Var g, C=0','Var g, rho=1.29','General Case','location','southwest');
hold off