function [t,y] = example2(tfinal,feedback,delay)

global a b r m vs

a = 0.7;
b = 0.8;
r = 0.08;
m = feedback;

vs = fzero(@vsstate,0)

history = [vs+0.1 0];
options = ddeset('RelTol',1e-8,'AbsTol',1e-9);

sol = dde23(@ddes,delay,history,[0,tfinal],options);

% Evaluate the solution for analysis and plotting.
t = linspace(0,tfinal,tfinal/0.5);
y = deval(sol,t)';

% Plot the solution.
figure(1)
plot(t,y(:,1),'k');
%axis([0 tfinal 0 6])

%------------------------------------------------------------------------
function dydt = ddes(t,y,ylag)

global a b r m vs

dydt = zeros(2,1);

v = y(1); 
w = y(2);

dydt(1) =  v - (1/3)*v^3 - w + m*(ylag(1) - vs);
dydt(2) =  r*(v + a - b*w);


%------------------------------------------------------------------------
function f = vsstate(v);

global a b

f = 1/3*v^3 + (1/b-1)*v + a/b;
