%
% REG with modifications by PHH 
% (note that must have lambda < mu -- see definitions)
%
% input epsilon (e.g. try 1, 0.2, 0.1, ...)
%
function MichMent
clear;
close;
time_end=5;
lambda=0.5;
mu=1.0;
%epsilon=0.01;
epsilon=input('epsilon=')
u=0.01:0.01:1; 
tau=(1-u-mu*log(u))/lambda;

options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4]);
[T,Y] = ode45(@mmeoms,[0 time_end],[1.0 0.0],options);
plot(T,Y(:,1),'-b',T,Y(:,2),'-r','LineWidth',3);
axis([0 time_end 0.0 1.0]);
hold on;
pause;
plot(T,Y(:,1)./(mu+Y(:,1)),'--r','LineWidth',3);
pause;
plot(tau,u,'--b','LineWidth',3);
pause;
hold off;
close;

function dy=mmeoms(t,y)
%  returns the coupled odes
    dy=zeros(2,1);
    dy(1)=-y(1)+(y(1)+mu-lambda)*y(2);
    dy(2)=(y(1)-(y(1)+mu)*y(2))/epsilon;
end
end







