function [t,altitude, velocity] = ode_parachute_jump_model()
	%% Parachute jumper model.
	% We assume that there is no x-component in the trajectory, 
    % i.e. the motion is only in the vertical direction.
    % 
	% We take in account the drag force which is directed against the motion
	% Fd =  -sign(v) * 1/2 * rho * v^2 * Cd * A
    % where
	% A   - the cross-sectional area of the object,
	% Cd  - the unit less drag coefficient (<=1) depends on the object shape,
	% rho - the density of the fluid,
	% v   - the velocity of the object

	% There is also the gravitational force m*g.
    % the resulting y acceleration is
	% y'' = -g  + Fd/m


	g=9.84;  % the acceleration due to gravity m/s^2
	rho=1.2; % kg/m^3 density of air
	Cd=.5;   % quite arbitrary choice
	m=100;   % the total mass of the jumper and the parachute in kg

	% transform system to the canonical form
	% t -> x
	% y -> y1
	% v = y' -> y2
	% thus
	% f1= y2
	% f2= -m*g  - sign(v)*|Fd| = -m*g  - sign(y2)*|Fd|

	function fvec = jumper_forces(x,y)
		fvec(1)=y(2);
		
		% We will use the time depended drag force.
		% For this problem it means that the parachute is closed during time Tc
		Tc=20; % seconds
		% The 'Aclosed' is  the area of the jumper + the closed parachute
		Aclosed=1; % m^2
		% after time Tc, the parachute opens and the combined area becomes much larger
		Aopen=60;  % m^2

		if (x<Tc) 
			% the parachute is closed
			A=Aclosed;
		else
			A=Aopen;
		end

		v=y(2); %velocity 

		Fd=1/2 * rho * v^2 * Cd * A; % the magnitude of the drag force
		fvec(2)= -g -sign(v)*Fd/m;

		% To make matlab happy, we need to return a column vector.
		% So, we transpose (note the dot in .')
		fvec=fvec.'; 
	end

	tspan=[0, 30]; % time interval of interest

	y0(1)=2000;    % initial altitude of the jumper
	y0(2)=0;       % initial velocity of the jumper

	% We are using MATLAB's ODE solver
	[t,y] = ode45(@jumper_forces, tspan, y0); 

    % Assigning the human readable variables names
	altitude=y(:,1); 
	velocity=y(:,2); 

	% The analytical no drag solution to provide a reference trajectory
	altitude_nodrag=y0(1)-g/2*t.^2;
	velocity_nodrag=y0(2)-g*t;

	figure(1);
	ax(1)=subplot(2,1,1);
	plot(t,altitude, 'r-', t(altitude_nodrag>0), altitude_nodrag(altitude_nodrag>0), 'b-');
	legend('with drag', 'no drag');
	xlabel('time');
	ylabel('altitude');

	ax(2)=subplot(2,1,2);
	plot(t,velocity, 'r-', t(altitude_nodrag>0), velocity_nodrag(altitude_nodrag>0), 'b-');
	legend('with drag', 'no drag');
	xlabel('time');
	ylabel('velocity');

	linkaxes(ax,'x'); % very handy for related subplots

end	
