% to use this, in the matlab command window, type
% odesolveresample;
function [v,w]=odesolverexample;
%  
% Solves system
%
% v'(t)= v*(a-v)*(v-1) - w + I(t);
% w'(t)= b*v - g*w;
%
% where a,b,g are constants, 0 < a < 1, and I(t) is a specified function.
%
% set constants
a = 0.7;
b = 0.01;
g = 0.03;
%
% specify initial conditions for v and w
%
v0 = 0;
w0 = 0.;
y0 = [v0 w0];
%
% define stimulus current
%
% Try I0 = 0.4, tstart = 0, tend = 1000
I0       = input(' Stimulus intensity (10) ')
tstart   = input(' Stimulus start time (0) ')
tend     = input(' Stimulus end time (10) ')
%
% define final time and frequency of outputs
%
% Try tfinal = 1500
tfinal   = input(' Final time (40) ')
% Try tstep = 0.5
tstep    = input(' Time step (0.1) ')
tspan    = [0:tstep:tfinal];
%
% solve ODEs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% call ode solver  -- stiff solver ode23s
%
%%%%%%%%%%%%%%%%%%
   [t,y] = ode23s(@odesys,tspan,y0);
%%%%%%%%%%%%%%%%%%
   tmax = max(t);
%
% plots of solutions
%
v = y(:,1);
w = y(:,2);
% 
figure
%
% plot v,w vs t
%
subplot(5,1,1:3)
   plot(t,v(:),'b-','LineWidth',2)
   ymin = min(v(:));
   ymax = max(v(:));
   ymin = ymin-0.5;
   ymax = ymax+0.5;
   axis([0 tmax ymin ymax])
   hold on
   plot(t,w(:),'r-','LineWidth',2)   
   xlabel('t','Fontsize',18)
   ylabel('v,w','Fontsize',18)
   legend('v','w')
   name1 = ' v and w vs t ----';
   name2 = strcat(' a=',num2str(a),',');
   name3 = strcat(' b=',num2str(b),',');
   name4 = strcat(' g=',num2str(g));
   name = strcat(name1,name2,name3,name4)
   title(name,'Fontsize',18)
% plot stimulus vs t   
jmax = length(t);
for j=1:jmax,
    time = t(j);
    stim(j) = getstim(time);
end    
subplot(5,1,5)
   numstep=length(tspan);
   iapp = zeros(numstep,1);
   plot(t,stim(:),'b-','Linewidth',2)
   ymin = min(stim(:)) - 0.1;
   ymax = max(stim(:)) + 0.1;
   axis([0 tmax ymin ymax])   
   xlabel('t','Fontsize',18)
   ylabel('I','Fontsize',18)
   title('I','Fontsize',18)   
   hold off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% begin nested function odesys
function dydt = odesys(t,y);
dydt = zeros(2,1);
v = y(1);
w = y(2);
I = getstim(t);
dydt(1) = v*(a-v)*(v-1) - w + I;
dydt(2) = b*v - g*w;
end  % end nested function odesys
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% begin nested function getstim
function I = getstim(t);
if t < tstart
   I = 0;
elseif tstart <= t & t <= tend
   I = I0;
else
   I = 0;
end
end  % end nested function getstim
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end  % end function odesolverexample