clear all

a= 3;   % prirodni prirastaj plena
b=1.1;  % stopa umanjenja plena usled grabljivaca
c=.5;    % natalitet grabljivaca
P0=1;
G0=0;
tkraj=5;

%P0=m/c; % pocetna populacija plena (stac.)
%G0=a/b; % pocetna populacija grabljivaca (stac.)


nizG = [0, 0.2, 0.5, 1, 2, 4, 6];
nizP = [0.1, 0.2, 0.5, 1, 2, 4, 6];

MaxP = zeros(length(nizG),length(nizP));
Steta = MaxP;

f=@(t,x) [a*x(1)-b*x(1)*x(2); c*x(1)];
for ii=1:length(nizG)
    G0 = nizG(ii);
    for jj=1:length(nizP)
        P0 = nizP(jj);
        [t,x]=ode45(f, [0,tkraj], [P0,G0]);
        MaxP(ii,jj) = max(x(:,1));
        Steta(ii,jj) = sum(x(:,1))/tkraj;
    end
end

figure(1)
surf( nizP, nizG, MaxP)
set(gca, 'YDir','reverse')
xlabel('P0 (virusi)')
ylabel('G0 (antitela)')
title('Max P (virusi)')
figure(2)
surf(nizP, nizG, Steta)
set(gca, 'YDir','reverse')
xlabel('P0 (virusi)')
ylabel('G0 (antitela)')
title('Steta po organizam (integral od P)')

