clear all
%Koliko stete virus napravi coveku?

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];

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

for i=1:length(nizG)
    G0 = nizG(i);
    for j=1:length(nizP)
        P0 = nizP(j);
        [t,x]=ode45(f, [0,tkraj], [P0,G0]);
        % Max kolicina virusa
        MaxP(i,j) = max(x(:,1));
        
        % Koliko je virus "izjedao" coveka (povrsina ispod krive)
        % Racunamo integral: I= (b-a)/2*(f(a)+f(b)) 
        % Tablica nije ekvidistantna pa za svaki interval odredjujemo h
        h=diff(t); %koraci za numericko integraljenje t_{i+1}-t_{i}
        zagrada=x(1:end-1,1)+x(2:end,1); %f_{i}+f_{i+1}
        Steta(i,j) = sum(h/2.*zagrada); %osnovna trapezna kvadraturna formula 
    end
end

% Maksimalne kolicine virusa za razlicite P0 i G0
figure(1)
surf( nizP, nizG, MaxP)
%hold on
%smrt=7;
%surf( nizP, nizG, repmat(smrt,length(nizP)))
set(gca, 'YDir','reverse')
xlabel('P0 (virusi)')
ylabel('G0 (antitela)')
title('Max P (virusi)')

% Koliko je covek ostecen kolicinski (povrsina ispod krive)
figure(2)
surf(nizP, nizG, Steta)
set(gca, 'YDir','reverse')
xlabel('P0 (virusi)')
ylabel('G0 (antitela)')
title('Steta po organizam (integral od P)')

