function NP
clear;
close;
time_end=200;
alpha=1.8;
beta=input('beta=')
delta=0.3;
uzero=0.5*(1-delta-alpha+sqrt((1-delta-alpha)^2+4*delta))
beta-(alpha*uzero^2/(uzero+delta)^2-uzero)
U=0:0.01:1.2;
V=(1-U).*(U+delta)/alpha;

options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4]);
[T,Y] = ode45(@lveoms,[0 time_end],[0.2 0.2],options);
plot(T,Y(:,1),'-b',T,Y(:,2),'-r','LineWidth',3);
axis([0 time_end 0.1 0.35]);
pause;
close;
[T,Y] = ode45(@lveoms,[0 time_end],[0.1 0.1],options);
plot(U,U,'--',U,V,'--','LineWidth',3);
axis([0 0.75 0 0.4]);
hold on;
pause;
plot(Y(:,1),Y(:,2),'-b','LineWidth',3);
pause;
[T,Y] = ode45(@lveoms,[0 time_end],[0.15 0.15],options);
plot(Y(:,1),Y(:,2),'-g','LineWidth',3);
pause;
[T,Y] = ode45(@lveoms,[0 time_end],[0.3 0.3],options);
plot(Y(:,1),Y(:,2),'-r','LineWidth',3);
pause;
hold off;
close;

function dy=lveoms(t,y)
%  returns the coupled odes
    dy=zeros(2,1);
    dy(1)=y(1)*(1-y(1))-alpha*y(1)*y(2)/(y(1)+delta);
    dy(2)=beta*y(2)*(1-y(2)/y(1));
end
end







 