function dy=raindropf(t,y);
global alpha beta gamma delta g

% No drag
%dy=[g-3*gamma*y(2)^(alpha-1)*y(1)^(beta+1);gamma*y(2)^alpha*y(1)^beta];

% Drag
dy=[g-3*gamma*y(2)^(alpha-1)*y(1)^(beta+1)-delta*(y(1)./y(2)).^(1.5); ...
gamma*y(2)^alpha*y(1)^beta];
