ex11ch4.m

these codes are for solving OED with matlab

function sol = ex11ch4
global a b r m v0
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(@ddes,20,history,[0, 60]);
   figure
   plot(sol.x,sol.y)
   title(['Plant''s Neuron Interaction Model' ...
          ' with m = ',num2str(m)])
   xlabel('time t')
end
%=====================================
function vp = ddes(t,y,Z)
global a b r m v0
ylag = Z(:,1);
v = y(1);
w = y(2);
vp = [ v - v^3 / 3 + w  + m * (ylag(1) - v0)
       r * (v + a - b*w)                     ];