function sol = prob3
% prob3, prob3f of tutorial.

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

% Compute steady state solution [v0; w0] for m = 0:
v0 = 2.417960226013935;
w0 = (v0 + a) / b;
history = [0.4*v0 ; 1.8*w0];

for m = [+10 -10]
   sol = dde23(@prob3f,20,history,[0, 60],[],v0,a,b,r,m);
   figure
   plot(sol.x,sol.y)
   title(['Problem 3. Plant''s Neuron Interaction Model' ...
          ' with m = ',num2str(m),'.'])
   xlabel('time t')
end

%=============================================================
function yp = prob3f(t,y,Z,v0,a,b,r,m)
ylag = Z(:,1);
v = y(1);
w = y(2);
yp = zeros(2,1);
yp(1) = v - v^3 / 3 + w  + m * (ylag(1) - v0);
yp(2) = r * (v + a - b*w);