function fall(gamma,u0,h0,N,T)
% integrates the free-fall equations via 1st order Euler's method
% dx2/dt = gee - gamma*x2; dx1/dt = x2;     x1: height (down is +)
% to run, start matlab, then type, e.g.:    x2: speed  ( "" "" "")
%               fall(.3,-10,0,1000,100)                 |
% gamma: air resistance         (sec^-1)                |
% u0   : initial velocity       (m/sec)         --------|----h=0
% h0   : initial height         (m)                     |
% T    : final time             (sec)                   | h,u
% N    : number of steps                                V
% gee  :acceleration of gravity (m/sec^2) 
eps = T/N;
x(1,1)=h0;x(1,2) = u0;gee=9.81;t(1) = 0;
for i = 1:N
    t(i+1)   = t(i)   + eps;
    x(i+1,2) = x(i,2) + eps*(gee - gamma*x(i,2));
    x(i+1,1) = x(i,1) + eps*x(i,2);
end
figure;
subplot(2,1,1)
axis([-2, 2 , -2, 2])
plot(t(:),-x(:,1),t(1:N/10:N),-x(1:N/10:N,1),'ro')
subplot(2,1,2)
axis([-2, 2 , -2, 2])
plot(t(:),-x(:,2),t(1:N/10:N),-x(1:N/10:N,2),'bo')
title(sprintf('phase plane for Free-Fall with eps = %12.5f',eps));

