function D = harmonic(position,t)
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % This function implements the differential equation %
   %                                                    %
   %              d/dt [x,y] = [-y,x].                  %
   %                                                    %
   % The analytic solutions are circles n the xy-plane: %
   %                                                    %
   %     [x,y] = A [cos t, sin t] + B [sin t, -cost].   %
   %                                                    %
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      
   x = position(1);
   y = position(2);
   D = [-y,x];
