function out = ML2
out{1} = @init;
out{2} = @fun_eval;
out{3} = [];
out{4} = [];
out{5} = [];
out{6} = [];
out{7} = [];
out{8} = [];
out{9} = [];

% --------------------------------------------------------------------------
function dydt = fun_eval(t,kmrgd,par_phi1,par_gCa1,par_V31,par_V41,par_Iapp1,par_phi2,par_gCa2,par_V32,par_V42,par_Iapp2,par_epsilon1,par_epsilon2)
CM=5;
gK=8;
gL=2;
VCa=120;
VK=-80;
VL=-60;
W1=-1.2;
W2=18;
minf1=(1+tanh((kmrgd(1)-W1)/W2))/2;
taun1=1/cosh((kmrgd(1)-par_V31)/(2*par_V41));
ninf1=(1+tanh((kmrgd(1)-par_V31)/par_V41))/2;
minf2=(1+tanh((kmrgd(3)-W1)/W2))/2;
taun2=1/cosh((kmrgd(3)-par_V32)/(2*par_V41));
ninf2=(1+tanh((kmrgd(3)-par_V32)/par_V41))/2;
dydt=[(-par_gCa1*minf1*(kmrgd(1)-VCa)-gK*kmrgd(2)*(kmrgd(1)-VK)-gL*(kmrgd(1)-VL)+par_Iapp1+par_epsilon1*(kmrgd(1)-kmrgd(3)))/CM;
par_phi1*(ninf1-kmrgd(2))/taun1;
(-par_gCa2*minf2*(kmrgd(3)-VCa)-gK*kmrgd(4)*(kmrgd(3)-VK)-gL*(kmrgd(3)-VL)+par_Iapp2+par_epsilon2*(kmrgd(3)-kmrgd(1)))/CM;
par_phi2*(ninf2-kmrgd(4))/taun2;];

% --------------------------------------------------------------------------
function [tspan,y0,options] = init
handles = feval(ML2);
y0=[0,0,0,0];
options = odeset('Jacobian',[],'JacobianP',[],'Hessians',[],'HessiansP',[]);
tspan = [0 10];

% --------------------------------------------------------------------------
function jac = jacobian(t,kmrgd,par_phi1,par_gCa1,par_V31,par_V41,par_Iapp1,par_phi2,par_gCa2,par_V32,par_V42,par_Iapp2,par_epsilon1,par_epsilon2)
% --------------------------------------------------------------------------
function jacp = jacobianp(t,kmrgd,par_phi1,par_gCa1,par_V31,par_V41,par_Iapp1,par_phi2,par_gCa2,par_V32,par_V42,par_Iapp2,par_epsilon1,par_epsilon2)
% --------------------------------------------------------------------------
function hess = hessians(t,kmrgd,par_phi1,par_gCa1,par_V31,par_V41,par_Iapp1,par_phi2,par_gCa2,par_V32,par_V42,par_Iapp2,par_epsilon1,par_epsilon2)
% --------------------------------------------------------------------------
function hessp = hessiansp(t,kmrgd,par_phi1,par_gCa1,par_V31,par_V41,par_Iapp1,par_phi2,par_gCa2,par_V32,par_V42,par_Iapp2,par_epsilon1,par_epsilon2)
%---------------------------------------------------------------------------
function tens3  = der3(t,kmrgd,par_phi1,par_gCa1,par_V31,par_V41,par_Iapp1,par_phi2,par_gCa2,par_V32,par_V42,par_Iapp2,par_epsilon1,par_epsilon2)
%---------------------------------------------------------------------------
function tens4  = der4(t,kmrgd,par_phi1,par_gCa1,par_V31,par_V41,par_Iapp1,par_phi2,par_gCa2,par_V32,par_V42,par_Iapp2,par_epsilon1,par_epsilon2)
%---------------------------------------------------------------------------
function tens5  = der5(t,kmrgd,par_phi1,par_gCa1,par_V31,par_V41,par_Iapp1,par_phi2,par_gCa2,par_V32,par_V42,par_Iapp2,par_epsilon1,par_epsilon2)
