function [t, rabbits, wolves] = ode_predator_prey_model();
	%% Predator prey model (Lotka-Volterra equation)
	% solving: 
	% dx/dt = x*(a - b*y)
	% dy/dt = y*(c*x - d)
	% x - prey (rabbits)  population
	% y - predator (wolves)  population
	% a - preys growth rate (we presume that they reproduce faster then die)
	% b - rate with which predators consumes preys
	% c - rate with which predators repopulate (proportional to catch size)
	% generally b ~= c
	% d - rate of death for predators 

	% first we transform to canonical form note that we have actually 2 1st order ODEs
	% t -> x
	% x -> y1
	% y -> y2 

	% so in matrix form it looks like
	% f1 =  y1' =  y1*(a - b*y2)
	% f2 =  y2' =  y2*(c*y1 - d)

	function fvec = predator_prey_forces(x,y)
		fvec(1)  =  y(1)*(a - b*y(2));
		fvec(2)  =  y(2)*(c*y(1) - d);
		% matlab compatibility issue we need to return a column vector
		% so we transpose (note the dot in .')
		fvec=fvec.';

	end

	% constants assignment
	a=  2; 	% rabbits repopulate quite fast
	b=1.0; 	% one wolf can kill a lot of rabbits
	c=.01;	% you need a lot of rabbits to give a live to a wolf
	d=.4; 	% wolves death rate

	% starting condition
	y0(1) =  115;  % starting size of the rabbits population
	y0(2) =    3;  % starting size of the wolves population

	tspan=[0, 100]; % remember 'x' has meaning of time

	% we will use matlab ODE solver
	options = odeset('RelTol',1e-8);
	[t, y]  = ode45(@predator_prey_forces, tspan, y0);

	rabbits=y(:,1);  
	wolves=y(:,2);

	figure(1);
	subplot(2,1,1);
	plot(t,rabbits, 'r-', t, wolves, 'b-');
	legend('rabbits', 'wolves');
	xlabel('time');
	ylabel('population size');

	subplot(2,1,2);
	plot(rabbits, wolves)
	xlabel('rabbits');
	ylabel('wolves');

end
