clear all

a= 3;   % prirodni prirastaj plena X
b=0.6;  % stopa umanjenja plena X usled grabljivaca Y
c=0.5;    % natalitet grabljivaca Y
d=0.3;  % mortalitet od Y
e= 0.5;    % stopa umanjenja Y usled grabljivaca Z
f= 0.2;    % natalitet grabljivaca Z
g= 1;    % mortalitet od Z

X0=50;
Y0=10;
Z0=5;

% Stacionarno:
%  g=a/3;  % Y=a/b=g/f - proizvoljno npr. da je odnos 1/3
%  f=b/3;
%  Y0=a/b;
%  d=c*X0-e*Z0;

f=@(t,x) [a*x(1)-b*x(1)*x(2); c*x(1)*x(2) - d*x(2)-e*x(2)*x(3); f*x(2)*x(3)-g*x(3)];

[t,x]=ode45(f, [0,5], [X0,Y0,Z0]);

figure(1)
%plot(t,x(:,1),'+',t,x(:,2),'*',t,x(:,3),'o');
plot(t,x(:,1),t,x(:,2),t,x(:,3));
xlabel('vreme')
ylabel('populacija')
legend('X','Y','Z')

figure(2)
plot3(x(:,1),x(:,2),x(:,3))
xlabel('X')
ylabel('Y')
zlabel('Z')
