global  ro_keha ro_ohk_0 g_rask csm mugrt
r_keha=0.1;
c_aero=0.47;
g_rask=9.814;
m_keha=15.;
p_ohk=101000.;
mu_ohk=29.;
r_gaas=8314.;
t_ohk=300.;

s_keha=pi*r_keha^2;
v_keha=4.*pi*r_keha^3/3.;
ro_ohk_0=p_ohk*mu_ohk/(r_gaas*t_ohk);
ro_keha=m_keha/v_keha;
mugrt=mu_ohk*g_rask/(r_gaas*t_ohk);
csm=0.5*c_aero*s_keha/m_keha;

y0=[10.;0.];
options=odeset('RelTol',1.e-6,'AbsTol',1.e-6,'MaxStep',1,'InitialStep',0.01,'Event',@mycheck) ;
[t,y,te,ye,ie]=ode45(@dim1_rof,[0 5],y0,options);

subplot (2,1,1)
plot(t,y(:,1),'*')
grid on
subplot(2,1,2)
plot(t,y(:,2))
grid on

function dy = dim1_rof(t,y)
global ro_keha ro_ohk_0 g_rask csm mugrt;
dy=zeros(2,1);
ro_ohk=ro_ohk_0*exp(-mugrt*y(1));
param1=g_rask*(ro_ohk/ro_keha-1.);
param2=csm*ro_ohk;
dy(1)=y(2);
dy(2)=param1-param2*abs(y(2))*y(2);
end

