t=0:0.05:3; %s
h0=1.5; %m
T0=298.15

function dhdt=dhdt(h)
   %Cantidades constantes:
   A=1; %m^2
   Cv1=0.02;
   Cv2=0.02;
   rho=1000;
   g=9.8;
   R=8.314;
   w=28.85;
   cv=722;


   %Parametros

   P1=150000;
   P3=100000;
   P0=100000;
   ng=150;
   M=ng*w*(1/1000);
   Vt=5;


   %Ecuaciones

   Vg=Vt-A*h(1);
   P0=ng*R*h(2)/Vg;
   P2=rho*g*h(1)+P0;
   F1=Cv1*sqrt(P1-P2);
   F2=Cv2*sqrt(P2-P3);
   dhdt(1)=(F1-F2)/A;
   dhdt(2)=(P0/(M*cv))*(F1-F2);

   endfunction


h=lsode("dhdt",[1.5,T0],t);