function [out1] = PROJECTILEfunction(t,y,flag,P)
% ----------------------------------------------------------------------
% o 2-D equations for projectile motion (see ex.4.3.2 from Fowles & Cassidy 2005)
% o x is the horizontal position, z the vertical position
% [see XX.m for further notes]
%   y(1) ... horiz. position x
%   y(2) ... horiz. velocity dx/dt
%   y(3) ... vert. position z
%   y(4) ... vert. velocity dz/dt

out1(1)= y(2);  % --> integrates to x(t)
out1(2)= -P.gamma* sqrt(y(2)^2 + y(4)^2)*y(2);  % --> integrates to dx/dt(t)
out1(3)= y(4);  % --> integrates to z(t)
out1(4)= -P.gamma* sqrt(y(2)^2 + y(4)^2)*y(4)- P.g;  % --> integrates to dz/dt(t)
out1= out1';    % wants output as a column vector