function Ydot = system(r,Y,omega)
    % M-file function for solution of the system of two differential equations
    Ydot = [ Y(2) ; -Y(2)/r - omega^2*Y(1) ];