function out = actinwavesys4again
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,k0,gam,k,del,s1,s2,r0,kn,ks,Tn,tau)
dydt=[(Tn-kmrgd(2))*(k0+gam*kmrgd(1)^3/(k^3+kmrgd(1)^3))-del*(s1+s2*kmrgd(3)/(r0+kmrgd(3)))*kmrgd(1);
(Tn-kmrgd(2))*(k0+gam*kmrgd(2)^3/(k^3+kmrgd(2)^3))-del*(s1+s2*kmrgd(4)/(r0+kmrgd(4)))*kmrgd(2);
tau*(kn*kmrgd(1)-ks*kmrgd(3));
tau*(kn*kmrgd(2)-ks*kmrgd(4));];

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

% --------------------------------------------------------------------------
function jac = jacobian(t,kmrgd,k0,gam,k,del,s1,s2,r0,kn,ks,Tn,tau)
% --------------------------------------------------------------------------
function jacp = jacobianp(t,kmrgd,k0,gam,k,del,s1,s2,r0,kn,ks,Tn,tau)
% --------------------------------------------------------------------------
function hess = hessians(t,kmrgd,k0,gam,k,del,s1,s2,r0,kn,ks,Tn,tau)
% --------------------------------------------------------------------------
function hessp = hessiansp(t,kmrgd,k0,gam,k,del,s1,s2,r0,kn,ks,Tn,tau)
%---------------------------------------------------------------------------
function tens3  = der3(t,kmrgd,k0,gam,k,del,s1,s2,r0,kn,ks,Tn,tau)
%---------------------------------------------------------------------------
function tens4  = der4(t,kmrgd,k0,gam,k,del,s1,s2,r0,kn,ks,Tn,tau)
%---------------------------------------------------------------------------
function tens5  = der5(t,kmrgd,k0,gam,k,del,s1,s2,r0,kn,ks,Tn,tau)
