function x_dot=rocket(t,x)
global a_o t_b c m_dot M_o M_f g_o a_f M_s gamma

v=x(1);
y=x(2);
rho=1.2*exp(-y/7000);
M=v/a_o;

if M <=1
    c_d=0.08*(1+exp(-(3*(1-M))^2));
else
    c_d=0.08*(1+exp(-1*(M-1)^2));
end

if t <=t_b
    T=m_dot*c;
    m=M_o-(m_dot*t);
else
    T=0;
    m=M_f;
end

D=c_d*.5*rho*(v^2)*a_f;

x_dot(1,1)=(T/m)-(D/m)-g_o;
x_dot(2,1)=v;