function LotkaVolterra()
%LOTKAVOLTERRA Summary of this function goes here
%   Detailed explanation goes here

	K = 2;
	a = 2;
	L = 3;
	b = 2;
	xmax = 3;
	ymax = 3;
	
	DF = @(x,y) [K - 2*x - a*y, -a*x;
		-b*y, L-b*x-2*y];
	
	nt = 1000;
	T = 20;
	t = linspace(0, T, nt);
	
	function [x,y] = simulateLV(t, x0, y0)
		x = zeros(size(t));
		y = zeros(size(t));
		x(1) = x0;
		y(1) = y0;
		dt = t(2)-t(1);
		nt = length(t);
		for n = 1:nt-1
			x(n+1) = x(n) + dt*x(n)*(K-x(n)-a*y(n));
			y(n+1) = y(n) + dt*y(n)*(L-b*x(n)-y(n));
		end
	end


	
	figure;
	plot([0, L/b], [L, 0], 'b', 'linewidth', 2);
	hold on;
	plot([0, xmax], [0, 0], 'b', 'linewidth', 2);
	plot([0, 0], [0, ymax], 'r', 'linewidth', 2);
	plot([0, K], [K/a, 0], 'r', 'linewidth', 2);
	
	u0 = [0 3
		1 3
		2 3
		3 3
		3 2
		3 0.5
		3 0
		1e-5 1e-9
		1e-5 1e-8
		1e-5 1e-7
		1e-5 1e-6
		1e-5 1e-5
		1e-6 1e-5];
	for i = 1 : size(u0,1)
		[x,y] = simulateLV(t, u0(i,1), u0(i,2));
		plot(x, y, 'k', 'linewidth', 2);
	end
	
	% The fixed point with xfp>0, yfp>0
	xfp = (K-a*L)/(1-a*b);
	yfp = (L-b*K)/(1-a*b);
	% Compute the eigenvectors of DF(xfp, yxp). The stable/unstable
	% subspaces will be parallel to rs and ru, respectively.
	[R,D] = eig(DF(xfp,yfp));
	rs = R(:,1);
	ru = R(:,2);
	[WsxMinus, WsyMinus]	= simulateLV(-t(1:260), xfp-rs(1)*1e-5, yfp-rs(2)*1e-5);
	[WsxPlus, WsyPlus]		= simulateLV(-t,		xfp+rs(1)*1e-5, yfp+rs(2)*1e-5);
	plot(WsxMinus, WsyMinus, 'g', 'linewidth', 2);
	plot(WsxPlus, WsyPlus, 'g', 'linewidth', 2);
	[WuxMinus, WuyMinus]	= simulateLV(t, xfp-ru(1)*1e-5, yfp-ru(2)*1e-5);
	[WuxPlus, WuyPlus]		= simulateLV(t, xfp+ru(1)*1e-5, yfp+ru(2)*1e-5);
	plot(WuxMinus, WuyMinus, 'm', 'linewidth', 2);
	plot(WuxPlus, WuyPlus, 'm', 'linewidth', 2);
	
	axis([0 xmax 0 ymax]);
end