% -- This Function outputs the derivatives of x, y and q --
% -- It is used in the main Program to solde the ODE --

function ds=derive(t, r, V) % r = [x y q]

rho = 0.5;
drho = 0.5;
                % r = c1 * v
q0 = 1;
k = 1850;

% Finding the Constants:
dx0 = 0.02;
Vs = 8.6;
ys = 0.04;
c1 = sqrt(drho / rho / ys);
qs = q0 / sqrt(dx0^2 + ys^2) / k;
c2 = c1 * c1 * dx0 * rho / qs / Vs;
% Derivations:

ds = r;

q = r(3);
x = r(1);
y = r(2);

rr = sqrt(x*x + y*y);

ds(3) = q0 - q * rr * k;
ds(1) = c2 * V * q / rho - c1 * c1 * x;
ds(2) = drho / rho - c1 * c1 * y * rr / 0.04;