function [t,altitude, velocity] = ode_parachute_jump_model()
	%% parachute jumper model
	% we solve some one jumping straight down (there is no x-component)
	% but we do not disregard the drag force which is against the motion direction
	% y'' = -g  + Fd/m

	% there is a drag force Fd in addition to the gravitational force m*g
	% Fd =  -sign(v) * 1/2 * rho * v^2 * Cd * A
	% A -  the cross-sectional area of the object, 
	% Cd - unit less drag coefficient (<=1) depends on the object shape, 
	% rho -  the density of the fluid,
	% v - speed of the object

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

	% 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 time depended drag force
		% for this problem it means that the parachute is closed during time Tc
		Tc=20; % seconds
		% and area of the jumper+parachute will be Ac
		Aclosed=1; % m^2
		% and then the parachute opens at area become 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;
		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

	% using matlab solver
	[t,y] = ode45(@jumper_forces, tspan, y0); 

	altitude=y(:,1);
	velocity=y(:,2); 

	% analytical no drag solution
	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');

end	
