% Here we look at a Morris-Lecar model 

clear

FunctionHandlesML = MorrisLecarModel;

% Numerically integrate the ODE system usig a stiff solver 
% If you look 2DMorrisLecarModel, you see that the 2nd part contains the function handle to the vector field
[t,y] = ode23s(FunctionHandlesML{2},[0 50],[0 0], [], 30, 10);

figure(1)
clf

% Plot the numerical solution 
plot(t, y);
title('Plot of y as a function of time');
xlabel('time'); ylabel('y(t)');

% Extract the initial equilibrium solution for solution continuation
x0 = y(end,:)';


% MATCONT settings
opt=contset;
opt=contset(opt,'Singularities',1);
opt=contset(opt,'MaxStepsize',10);
opt=contset(opt,'MaxNumpoints',200);

% initialization call 
[x1,v1]=init_EP_EP(@MorrisLecarModel,x0,[30;6],[1]);

% continuation
[x,v,s,h,f]=cont(@equilibrium,x1,[],opt);


figure(2)
clf

% Plot the solution 
cpl(x,v,s,[(length(x0)+1) 1]);
axis([0 400 -50 50])