function [t, x, y, vx, vy] = ode_projectile_with_air_drag_model()
	%% solving the equation of motions for a projectile with air drug included
	% r''= F = Fg+ Fd,  Fg is a gravity pull, Fd is the air drug force
	% this can be decomposed to
	% x'' = Fd_x
	% y'' = -g + Fd_y
	% for more details about air drug see ode_parachute_jump_model.m
	% Fd =  1/2 * rho * v^2 * Cd * A
	% thus
	% Fd_x= - Fd * v_x/v  ; % vx/v takes care of the proper sign of the drag projection
	% Fd_y= - Fd * v_y/v  ; % vy/v takes care of the proper sign of the drag projection

	% at the first look it does not look like ODE but since x and y depends only on t 
	% it is actually a system of ODEs

	% transform system to the canonical form
	% t  -> x
	% x  -> y1
	% vx -> y2
	% y  -> y3
	% vy -> y4
	%
	% f1 -> y2
	% f2 -> Fd_x
	% f3 -> y4
	% f4 -> -g + Fd_y

	% some constants 
	rho=1.2; % kg/m^3
	Cd=.5;   % quite arbitrary choice
	m=0.01; % mass of the projectile  and the parachute in kg
	g=9.84; % acceleration due to gravity
	A=1e-4; % in m^2  , quite typical for a bullet

	function fvec = projectile_forces(x,y)
		vx=y(2);
		vy=y(4);
		v=sqrt(vx^2+vy^2);

		Fd=1/2 * rho * v^2 * Cd * A;

		fvec(1) = y(2);
		fvec(2) = -Fd*vx/v/m;
		fvec(3) = y(4);
		fvec(4) = -g -Fd*vy/v/m;

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

	tspan=[0, 80]; % time interval of interest
	theta=pi/4;    % shooting angle
	v0 = 700;      % initial projectile speed in m/s (quite typical for a bullet)
	y0(1)=0;   % initial x position
	y0(2)=v0*cos(theta); % vx
	y0(3)=0;   % initial y position
	y0(4)=v0*sin(theta); % vy

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

	x =  ysol(:,1);
	vx = ysol(:,2);
	y =  ysol(:,3);
	vy = ysol(:,4);

	x_analytical = y0(1) + y0(2)*t;
	y_analytical = y0(3) + y0(4)*t -g/2*t.^2;
	figure(1);
	subplot(2,1,1);
	plot(x,y, 'r-', x_analytical, y_analytical, 'b-');
	legend('with drag', 'no drag');
	xlabel('x');
	ylabel('y');
	title ('y vs x');

	v=sqrt(vx.^2+vy.^2);
	v_analytical= sqrt(y0(2).^2 + (y0(4) - g*t).^2) ;
	subplot(2,1,2);
	plot(x,v, 'r-', x_analytical, v_analytical, 'b-');
	legend('with drag', 'no drag');
	xlabel('x');
	ylabel('v');
	title ('speed vs x');
end
