function harvest()
clear all

K=350;  % maksimalni broj jedinki
N0=300; % pocetni broj jedinki u populaciji
r=0.09; % prirastaj r=p(0)
%H=1/4*r*K % kriticna zetva (harvest)
H=1/4*r*K-1 % stabilna
%H=6

disp('N1 i N2')
roots([1 -K H*K/r])

% Ako zelimo da prekinemo ode45 u trenutku kada se dostigne x=200 pravimo Event
%Opt    = odeset('Events', @myEvent);
%[t,x]=ode45(@(t,x) x*r*(1-x/K)-H,[0,100],N0,Opt);

[t,x]=ode45(@(t,x) x*r*(1-x/K)-H,[0,500],N0);

%t(end)
%x(end)
plot(t,x)

function [value,isterminal,direction] = myEvent(t,x)
value = x-200; % The value that we want to be zero
isterminal = 1;  % Halt integration 
direction = 0;   % The zero can be approached from either direction
