global  ro_keha ro_ohk_0 g_rask param1 param2
r_keha=0.5;
c_aero=0.47;
g_rask=9.814;
m_keha=15.;
p_ohk=102500.;
mu_ohk=29.;
r_gaas=8314.;
t_ohk=266.;

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;
param2=mu_ohk*g_rask/(r_gaas*t_ohk);
param1=0.5*c_aero*s_keha/m_keha;

y0=[0.;0.;10.;20.];
options=odeset('RelTol',1.e-8,'AbsTol',1.e-8,'MaxStep',0.001,'InitialStep',0.0011,'Event',@mycheck) ;
[t,y,te,ye,ie]=ode45(@dim2_rof,[0 5],y0,options);
subplot (2,1,1)
plot(y(:,1),y(:,2))
grid on
subplot(2,1,2)
plot(y(:,2),sqrt(y(:,3).^2+y(:,4).^2))

grid on

function dy = dim2_rof(t,y)
global ro_keha ro_ohk_0 g_rask param1 param2;
dy=zeros(4,1);
v_abs=sqrt(y(3)*y(3)+y(4)*y(4));
ro_ohk=ro_ohk_0*exp(-param2*y(2));
param3=g_rask*(ro_ohk/ro_keha-1.);
param4=param1*ro_ohk;
dy(1)=y(3);
dy(2)=y(4);
dy(3)=-param4*v_abs*y(3);
dy(4)=param3-param4*v_abs*y(4);
end

