h =  linspace(0,40,1000);
D = density2(h*1000);
plot(h,D,'k')
xlabel('h (km)')
ylabel('\rho (kg/m^3)')