function out = forcedML2
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_a,par_omega,par_epsilon,par_phi1,par_phi2,par_CM,par_gCa1,par_gCa2,par_V31,par_V32,par_V41,par_V42)
gK=8;
gL=2;
VCa=120;
VK=-80;
VL=-60;
W1=-1.2;
W2=18;
minf1=(1+tanh((kmrgd(3)-W1)/W2))/2;
taun1=1/cosh((kmrgd(3)-par_V31)/(2*par_V41));
ninf1=(1+tanh((kmrgd(3)-par_V31)/par_V41))/2;
minf2=(1+tanh((kmrgd(5)-W1)/W2))/2;
taun2=1/cosh((kmrgd(5)-par_V32)/(2*par_V42));
ninf2=(1+tanh((kmrgd(5)-par_V32)/par_V42))/2;
dydt=[par_a*kmrgd(1)-par_omega*kmrgd(2)-kmrgd(1)*(kmrgd(1)^2+kmrgd(2)^2);
par_omega*kmrgd(1)+par_a*kmrgd(2)-kmrgd(2)*(kmrgd(1)^2+kmrgd(2)^2);
(-par_gCa1*minf1*(kmrgd(3)-VCa)-gK*kmrgd(4)*(kmrgd(3)-VK)-gL*(kmrgd(3)-VL)+kmrgd(1)+par_epsilon*(kmrgd(3)-kmrgd(5)))/par_CM;
par_phi1*(ninf1-kmrgd(4))/taun1;
(-par_gCa2*minf2*(kmrgd(5)-VCa)-gK*kmrgd(6)*(kmrgd(5)-VK)-gL*(kmrgd(5)-VL)+kmrgd(1)+par_epsilon*(kmrgd(5)-kmrgd(3)))/par_CM;
par_phi2*(ninf2-kmrgd(6))/taun2;];

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

% --------------------------------------------------------------------------
function jac = jacobian(t,kmrgd,par_a,par_omega,par_epsilon,par_phi1,par_phi2,par_CM,par_gCa1,par_gCa2,par_V31,par_V32,par_V41,par_V42)
% --------------------------------------------------------------------------
function jacp = jacobianp(t,kmrgd,par_a,par_omega,par_epsilon,par_phi1,par_phi2,par_CM,par_gCa1,par_gCa2,par_V31,par_V32,par_V41,par_V42)
% --------------------------------------------------------------------------
function hess = hessians(t,kmrgd,par_a,par_omega,par_epsilon,par_phi1,par_phi2,par_CM,par_gCa1,par_gCa2,par_V31,par_V32,par_V41,par_V42)
% --------------------------------------------------------------------------
function hessp = hessiansp(t,kmrgd,par_a,par_omega,par_epsilon,par_phi1,par_phi2,par_CM,par_gCa1,par_gCa2,par_V31,par_V32,par_V41,par_V42)
%---------------------------------------------------------------------------
function tens3  = der3(t,kmrgd,par_a,par_omega,par_epsilon,par_phi1,par_phi2,par_CM,par_gCa1,par_gCa2,par_V31,par_V32,par_V41,par_V42)
%---------------------------------------------------------------------------
function tens4  = der4(t,kmrgd,par_a,par_omega,par_epsilon,par_phi1,par_phi2,par_CM,par_gCa1,par_gCa2,par_V31,par_V32,par_V41,par_V42)
%---------------------------------------------------------------------------
function tens5  = der5(t,kmrgd,par_a,par_omega,par_epsilon,par_phi1,par_phi2,par_CM,par_gCa1,par_gCa2,par_V31,par_V32,par_V41,par_V42)
