    function dydt=rhs(t,y)
    
       global A T Ku Kl m b Cd rho Ac ;

    
       
       x = A*cos(2*pi*t/T);
       x = x*(1-exp(-.01*t.*t.*t/T));
       
       %b
       %foo = .5*Cd*rho*Ac*abs(y(2))
       
       
       dydt_1 = y(2);
       %dydt_2 = (-(Ku+Kl)*y(1) - b*y(2) + Ku*x)/m;  %Linear Damping
       dydt_2 = (-(Ku+Kl)*y(1) - 0.5*Cd*rho*Ac*y(2)*abs(y(2)) + Ku*x)/m;  %Quadratic Damping
       
       dydt=[dydt_1; dydt_2];
    end
