function dy = ballf(t,y)
global g CD CL alpha
dy = zeros(4,1);    % a column vector
v = sqrt(y(3).^2+y(4).^2);   % speed v
dy(1) = y(3);
dy(2) = y(4);
dy(3) = -alpha*v.*(CD*y(3)+CL*y(4));
dy(4) = alpha*v.*(-CD*y(4)+CL*y(3))-g;