function out = forcedML
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_phi,par_CM,par_gCa,par_V3,par_V4)
gK=8;
gL=2;
VCa=120;
VK=-80;
VL=-60;
V1=-1.2;
V2=18;
minf=(1+tanh((kmrgd(3)-V1)/V2))/2;
taun=1/cosh((kmrgd(3)-par_V3)/(2*par_V4));
ninf=(1+tanh((kmrgd(3)-par_V3)/par_V4))/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_gCa*minf*(kmrgd(3)-VCa)-gK*kmrgd(4)*(kmrgd(3)-VK)-gL*(kmrgd(3)-VL)+kmrgd(1))/par_CM;
par_phi*(ninf-kmrgd(4))/taun;];

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

% --------------------------------------------------------------------------
function jac = jacobian(t,kmrgd,par_a,par_omega,par_phi,par_CM,par_gCa,par_V3,par_V4)
% --------------------------------------------------------------------------
function jacp = jacobianp(t,kmrgd,par_a,par_omega,par_phi,par_CM,par_gCa,par_V3,par_V4)
% --------------------------------------------------------------------------
function hess = hessians(t,kmrgd,par_a,par_omega,par_phi,par_CM,par_gCa,par_V3,par_V4)
% --------------------------------------------------------------------------
function hessp = hessiansp(t,kmrgd,par_a,par_omega,par_phi,par_CM,par_gCa,par_V3,par_V4)
%---------------------------------------------------------------------------
function tens3  = der3(t,kmrgd,par_a,par_omega,par_phi,par_CM,par_gCa,par_V3,par_V4)
%---------------------------------------------------------------------------
function tens4  = der4(t,kmrgd,par_a,par_omega,par_phi,par_CM,par_gCa,par_V3,par_V4)
%---------------------------------------------------------------------------
function tens5  = der5(t,kmrgd,par_a,par_omega,par_phi,par_CM,par_gCa,par_V3,par_V4)
