function out = PVIN_Cai
out{1} = @init;
out{2} = @fun_eval;
out{3} = @jacobian;
out{4} = @jacobianp;
out{5} = @hessians;
out{6} = @hessiansp;
out{7} = @der3;
out{8} = [];
out{9} = [];

% --------------------------------------------------------------------------
function dydt = fun_eval(t,kmrgd,par_Cai,par_Iapp)
gNa=300;
VNa=58;
gKv1=15;
gKv3=180;
VK=-80;
gCa=8;
VCa=68;
Vleak=-50;
gleak=8;
gSK=10;
Cm=30;
nk=5;
ksk=0.8;
pgamma=0.01;
mArea=3000;
Vm=-17.5;
Sm=-11.4;
Aah=0.0025;
Sah=10;
Vah=23;
Abh=0.094;
Sbh=-5.5;
Vbh=-31;
Aan1=0.002;
San1=-9;
Van1=-30;
Abn1=0.017;
Vbn1=-35;
Sbn1=5.9;
Aan3=1.98;
San3=-12.6;
Van3=96;
Abn3=0.34;
Vbn3=-36;
Sbn3=10.5;
Va=3;
Sa=-10.4;
Bt=0;
mmax=1/(1+exp((kmrgd(1)-Vm)/Sm));
ah=Aah/exp((kmrgd(1)-Vah)/Sah);
bh=Abh*(kmrgd(1)-Vbh)/(1-exp((kmrgd(1)-Vbh)/Sbh));
amax=1/(1+exp((kmrgd(1)-Va)/Sa));
an1=Aan1*(kmrgd(1)-Van1)/(1-exp((kmrgd(1)-Van1)/San1));
bn1=Abn1/exp((kmrgd(1)-Vbn1)/Sbn1);
an3=Aan3*(kmrgd(1)-Van3)/(1-exp((kmrgd(1)-Van3)/San3));
bn3=Abn3/exp((kmrgd(1)-Vbn3)/Sbn3);
k=par_Cai^nk/(ksk^nk+par_Cai^nk);
dydt=[(-gNa*mmax^3*kmrgd(2)*(kmrgd(1)-VNa)-gKv1*kmrgd(3)^4*(kmrgd(1)-VK)-gKv3*kmrgd(4)^2*(kmrgd(1)-VK)-gCa*amax^2*(kmrgd(1)-VCa)-gSK*k*(kmrgd(1)-VK)-gleak*(kmrgd(1)-Vleak)+par_Iapp)/Cm;
ah*(1-kmrgd(2))-bh*kmrgd(2);
an1*(1-kmrgd(3))-bn1*kmrgd(3);
an3*(1-kmrgd(4))-bn3*kmrgd(4);];

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

% --------------------------------------------------------------------------
function jac = jacobian(t,kmrgd,par_Cai,par_Iapp)
jac=[ - (10*kmrgd(2))/(exp(- (5*kmrgd(1))/57 - 175/114) + 1)^3 - 4/(15*(exp(15/52 - (5*kmrgd(1))/52) + 1)^2) - par_Cai^5/(3*(par_Cai^5 + 1024/3125)) - kmrgd(3)^4/2 - 6*kmrgd(4)^2 - (2*exp(15/52 - (5*kmrgd(1))/52)*(kmrgd(1) - 68))/(39*(exp(15/52 - (5*kmrgd(1))/52) + 1)^3) - (50*kmrgd(2)*exp(- (5*kmrgd(1))/57 - 175/114)*(kmrgd(1) - 58))/(19*(exp(- (5*kmrgd(1))/57 - 175/114) + 1)^4) - 4/15 , -(10*(kmrgd(1) - 58))/(exp(- (5*kmrgd(1))/57 - 175/114) + 1)^3 , -2*kmrgd(3)^3*(kmrgd(1) + 80) , -12*kmrgd(4)*(kmrgd(1) + 80) ; (47*kmrgd(2))/(500*(exp(- (2*kmrgd(1))/11 - 62/11) - 1)) + (exp(23/10 - kmrgd(1)/10)*(kmrgd(2) - 1))/4000 + (2*kmrgd(2)*exp(- (2*kmrgd(1))/11 - 62/11)*((47*kmrgd(1))/500 + 1457/500))/(11*(exp(- (2*kmrgd(1))/11 - 62/11) - 1)^2) , ((47*kmrgd(1))/500 + 1457/500)/(exp(- (2*kmrgd(1))/11 - 62/11) - 1) - exp(23/10 - kmrgd(1)/10)/400 , 0 , 0 ; (17*kmrgd(3)*exp(- (10*kmrgd(1))/59 - 350/59))/5900 + (kmrgd(3) - 1)/(500*(exp(- kmrgd(1)/9 - 10/3) - 1)) + (exp(- kmrgd(1)/9 - 10/3)*(kmrgd(1)/500 + 3/50)*(kmrgd(3) - 1))/(9*(exp(- kmrgd(1)/9 - 10/3) - 1)^2) , 0 , (kmrgd(1)/500 + 3/50)/(exp(- kmrgd(1)/9 - 10/3) - 1) - (17*exp(- (10*kmrgd(1))/59 - 350/59))/1000 , 0 ; (17*kmrgd(4)*exp(- (2*kmrgd(1))/21 - 24/7))/525 + (99*(kmrgd(4) - 1))/(50*(exp(160/21 - (5*kmrgd(1))/63) - 1)) + (5*exp(160/21 - (5*kmrgd(1))/63)*((99*kmrgd(1))/50 - 4752/25)*(kmrgd(4) - 1))/(63*(exp(160/21 - (5*kmrgd(1))/63) - 1)^2) , 0 , 0 , ((99*kmrgd(1))/50 - 4752/25)/(exp(160/21 - (5*kmrgd(1))/63) - 1) - (17*exp(- (2*kmrgd(1))/21 - 24/7))/50 ];
% --------------------------------------------------------------------------
function jacp = jacobianp(t,kmrgd,par_Cai,par_Iapp)
jacp=[ (5*par_Cai^9*(kmrgd(1) + 80))/(3*(par_Cai^5 + 1024/3125)^2) - (5*par_Cai^4*(kmrgd(1) + 80))/(3*(par_Cai^5 + 1024/3125)) , 1/30 ; 0 , 0 ; 0 , 0 ; 0 , 0 ];
% --------------------------------------------------------------------------
function hess = hessians(t,kmrgd,par_Cai,par_Iapp)
hess1=[ (5*exp(15/52 - (5*kmrgd(1))/52)*(kmrgd(1) - 68))/(1014*(exp(15/52 - (5*kmrgd(1))/52) + 1)^3) - (100*kmrgd(2)*exp(- (5*kmrgd(1))/57 - 175/114))/(19*(exp(- (5*kmrgd(1))/57 - 175/114) + 1)^4) - (5*exp(15/26 - (5*kmrgd(1))/26)*(kmrgd(1) - 68))/(338*(exp(15/52 - (5*kmrgd(1))/52) + 1)^4) - (4*exp(15/52 - (5*kmrgd(1))/52))/(39*(exp(15/52 - (5*kmrgd(1))/52) + 1)^3) - (1000*kmrgd(2)*exp(- (10*kmrgd(1))/57 - 175/57)*(kmrgd(1) - 58))/(1083*(exp(- (5*kmrgd(1))/57 - 175/114) + 1)^5) + (250*kmrgd(2)*exp(- (5*kmrgd(1))/57 - 175/114)*(kmrgd(1) - 58))/(1083*(exp(- (5*kmrgd(1))/57 - 175/114) + 1)^4) , - 10/(exp(- (5*kmrgd(1))/57 - 175/114) + 1)^3 - (50*exp(- (5*kmrgd(1))/57 - 175/114)*(kmrgd(1) - 58))/(19*(exp(- (5*kmrgd(1))/57 - 175/114) + 1)^4) , -2*kmrgd(3)^3 , -12*kmrgd(4) ; (47*kmrgd(2)*exp(- (2*kmrgd(1))/11 - 62/11))/(1375*(exp(- (2*kmrgd(1))/11 - 62/11) - 1)^2) - (exp(23/10 - kmrgd(1)/10)*(kmrgd(2) - 1))/40000 - (4*kmrgd(2)*exp(- (2*kmrgd(1))/11 - 62/11)*((47*kmrgd(1))/500 + 1457/500))/(121*(exp(- (2*kmrgd(1))/11 - 62/11) - 1)^2) + (8*kmrgd(2)*exp(- (4*kmrgd(1))/11 - 124/11)*((47*kmrgd(1))/500 + 1457/500))/(121*(exp(- (2*kmrgd(1))/11 - 62/11) - 1)^3) , exp(23/10 - kmrgd(1)/10)/4000 + 47/(500*(exp(- (2*kmrgd(1))/11 - 62/11) - 1)) + (2*exp(- (2*kmrgd(1))/11 - 62/11)*((47*kmrgd(1))/500 + 1457/500))/(11*(exp(- (2*kmrgd(1))/11 - 62/11) - 1)^2) , 0 , 0 ; (exp(- kmrgd(1)/9 - 10/3)*(kmrgd(3) - 1))/(2250*(exp(- kmrgd(1)/9 - 10/3) - 1)^2) - (17*kmrgd(3)*exp(- (10*kmrgd(1))/59 - 350/59))/34810 - (exp(- kmrgd(1)/9 - 10/3)*(kmrgd(1)/500 + 3/50)*(kmrgd(3) - 1))/(81*(exp(- kmrgd(1)/9 - 10/3) - 1)^2) + (2*exp(- (2*kmrgd(1))/9 - 20/3)*(kmrgd(1)/500 + 3/50)*(kmrgd(3) - 1))/(81*(exp(- kmrgd(1)/9 - 10/3) - 1)^3) , 0 , (17*exp(- (10*kmrgd(1))/59 - 350/59))/5900 + 1/(500*(exp(- kmrgd(1)/9 - 10/3) - 1)) + (exp(- kmrgd(1)/9 - 10/3)*(kmrgd(1)/500 + 3/50))/(9*(exp(- kmrgd(1)/9 - 10/3) - 1)^2) , 0 ; (11*exp(160/21 - (5*kmrgd(1))/63)*(kmrgd(4) - 1))/(35*(exp(160/21 - (5*kmrgd(1))/63) - 1)^2) - (34*kmrgd(4)*exp(- (2*kmrgd(1))/21 - 24/7))/11025 - (25*exp(160/21 - (5*kmrgd(1))/63)*((99*kmrgd(1))/50 - 4752/25)*(kmrgd(4) - 1))/(3969*(exp(160/21 - (5*kmrgd(1))/63) - 1)^2) + (50*exp(320/21 - (10*kmrgd(1))/63)*((99*kmrgd(1))/50 - 4752/25)*(kmrgd(4) - 1))/(3969*(exp(160/21 - (5*kmrgd(1))/63) - 1)^3) , 0 , 0 , (17*exp(- (2*kmrgd(1))/21 - 24/7))/525 + 99/(50*(exp(160/21 - (5*kmrgd(1))/63) - 1)) + (5*exp(160/21 - (5*kmrgd(1))/63)*((99*kmrgd(1))/50 - 4752/25))/(63*(exp(160/21 - (5*kmrgd(1))/63) - 1)^2) ];
hess2=[ - 10/(exp(- (5*kmrgd(1))/57 - 175/114) + 1)^3 - (50*exp(- (5*kmrgd(1))/57 - 175/114)*(kmrgd(1) - 58))/(19*(exp(- (5*kmrgd(1))/57 - 175/114) + 1)^4) , 0 , 0 , 0 ; exp(23/10 - kmrgd(1)/10)/4000 + 47/(500*(exp(- (2*kmrgd(1))/11 - 62/11) - 1)) + (2*exp(- (2*kmrgd(1))/11 - 62/11)*((47*kmrgd(1))/500 + 1457/500))/(11*(exp(- (2*kmrgd(1))/11 - 62/11) - 1)^2) , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ];
hess3=[ -2*kmrgd(3)^3 , 0 , -6*kmrgd(3)^2*(kmrgd(1) + 80) , 0 ; 0 , 0 , 0 , 0 ; (17*exp(- (10*kmrgd(1))/59 - 350/59))/5900 + 1/(500*(exp(- kmrgd(1)/9 - 10/3) - 1)) + (exp(- kmrgd(1)/9 - 10/3)*(kmrgd(1)/500 + 3/50))/(9*(exp(- kmrgd(1)/9 - 10/3) - 1)^2) , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ];
hess4=[ -12*kmrgd(4) , 0 , 0 , - 12*kmrgd(1) - 960 ; 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ; (17*exp(- (2*kmrgd(1))/21 - 24/7))/525 + 99/(50*(exp(160/21 - (5*kmrgd(1))/63) - 1)) + (5*exp(160/21 - (5*kmrgd(1))/63)*((99*kmrgd(1))/50 - 4752/25))/(63*(exp(160/21 - (5*kmrgd(1))/63) - 1)^2) , 0 , 0 , 0 ];
hess(:,:,1) =hess1;
hess(:,:,2) =hess2;
hess(:,:,3) =hess3;
hess(:,:,4) =hess4;
% --------------------------------------------------------------------------
function hessp = hessiansp(t,kmrgd,par_Cai,par_Iapp)
hessp1=[ (5*par_Cai^9)/(3*(par_Cai^5 + 1024/3125)^2) - (5*par_Cai^4)/(3*(par_Cai^5 + 1024/3125)) , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ];
hessp2=[ 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ];
hessp(:,:,1) =hessp1;
hessp(:,:,2) =hessp2;
%---------------------------------------------------------------------------
function tens3  = der3(t,kmrgd,par_Cai,par_Iapp)
tens31=[ (5*exp(15/52 - (5*kmrgd(1))/52))/(338*(exp(15/52 - (5*kmrgd(1))/52) + 1)^3) - (15*exp(15/26 - (5*kmrgd(1))/26))/(338*(exp(15/52 - (5*kmrgd(1))/52) + 1)^4) - (1000*kmrgd(2)*exp(- (10*kmrgd(1))/57 - 175/57))/(361*(exp(- (5*kmrgd(1))/57 - 175/114) + 1)^5) + (250*kmrgd(2)*exp(- (5*kmrgd(1))/57 - 175/114))/(361*(exp(- (5*kmrgd(1))/57 - 175/114) + 1)^4) + (75*exp(15/26 - (5*kmrgd(1))/26)*(kmrgd(1) - 68))/(17576*(exp(15/52 - (5*kmrgd(1))/52) + 1)^4) - (25*exp(15/52 - (5*kmrgd(1))/52)*(kmrgd(1) - 68))/(52728*(exp(15/52 - (5*kmrgd(1))/52) + 1)^3) + (5000*kmrgd(2)*exp(- (10*kmrgd(1))/57 - 175/57)*(kmrgd(1) - 58))/(20577*(exp(- (5*kmrgd(1))/57 - 175/114) + 1)^5) - (1250*kmrgd(2)*exp(- (5*kmrgd(1))/57 - 175/114)*(kmrgd(1) - 58))/(61731*(exp(- (5*kmrgd(1))/57 - 175/114) + 1)^4) - (25*exp(15/26 - (5*kmrgd(1))/26)*exp(15/52 - (5*kmrgd(1))/52)*(kmrgd(1) - 68))/(4394*(exp(15/52 - (5*kmrgd(1))/52) + 1)^5) - (25000*kmrgd(2)*exp(- (10*kmrgd(1))/57 - 175/57)*exp(- (5*kmrgd(1))/57 - 175/114)*(kmrgd(1) - 58))/(61731*(exp(- (5*kmrgd(1))/57 - 175/114) + 1)^6) , (250*exp(- (5*kmrgd(1))/57 - 175/114)*(kmrgd(1) - 58))/(1083*(exp(- (5*kmrgd(1))/57 - 175/114) + 1)^4) - (1000*exp(- (10*kmrgd(1))/57 - 175/57)*(kmrgd(1) - 58))/(1083*(exp(- (5*kmrgd(1))/57 - 175/114) + 1)^5) - (100*exp(- (5*kmrgd(1))/57 - 175/114))/(19*(exp(- (5*kmrgd(1))/57 - 175/114) + 1)^4) , 0 , 0 ; (exp(23/10 - kmrgd(1)/10)*(kmrgd(2) - 1))/400000 - (141*kmrgd(2)*exp(- (2*kmrgd(1))/11 - 62/11))/(15125*(exp(- (2*kmrgd(1))/11 - 62/11) - 1)^2) + (282*kmrgd(2)*exp(- (4*kmrgd(1))/11 - 124/11))/(15125*(exp(- (2*kmrgd(1))/11 - 62/11) - 1)^3) + (8*kmrgd(2)*exp(- (2*kmrgd(1))/11 - 62/11)*((47*kmrgd(1))/500 + 1457/500))/(1331*(exp(- (2*kmrgd(1))/11 - 62/11) - 1)^2) - (48*kmrgd(2)*exp(- (4*kmrgd(1))/11 - 124/11)*((47*kmrgd(1))/500 + 1457/500))/(1331*(exp(- (2*kmrgd(1))/11 - 62/11) - 1)^3) + (48*kmrgd(2)*exp(- (2*kmrgd(1))/11 - 62/11)*exp(- (4*kmrgd(1))/11 - 124/11)*((47*kmrgd(1))/500 + 1457/500))/(1331*(exp(- (2*kmrgd(1))/11 - 62/11) - 1)^4) , (47*exp(- (2*kmrgd(1))/11 - 62/11))/(1375*(exp(- (2*kmrgd(1))/11 - 62/11) - 1)^2) - exp(23/10 - kmrgd(1)/10)/40000 - (4*exp(- (2*kmrgd(1))/11 - 62/11)*((47*kmrgd(1))/500 + 1457/500))/(121*(exp(- (2*kmrgd(1))/11 - 62/11) - 1)^2) + (8*exp(- (4*kmrgd(1))/11 - 124/11)*((47*kmrgd(1))/500 + 1457/500))/(121*(exp(- (2*kmrgd(1))/11 - 62/11) - 1)^3) , 0 , 0 ; (17*kmrgd(3)*exp(- (10*kmrgd(1))/59 - 350/59))/205379 - (exp(- kmrgd(1)/9 - 10/3)*(kmrgd(3) - 1))/(13500*(exp(- kmrgd(1)/9 - 10/3) - 1)^2) + (exp(- (2*kmrgd(1))/9 - 20/3)*(kmrgd(3) - 1))/(6750*(exp(- kmrgd(1)/9 - 10/3) - 1)^3) + (exp(- kmrgd(1)/9 - 10/3)*(kmrgd(1)/500 + 3/50)*(kmrgd(3) - 1))/(729*(exp(- kmrgd(1)/9 - 10/3) - 1)^2) - (2*exp(- (2*kmrgd(1))/9 - 20/3)*(kmrgd(1)/500 + 3/50)*(kmrgd(3) - 1))/(243*(exp(- kmrgd(1)/9 - 10/3) - 1)^3) + (2*exp(- kmrgd(1)/9 - 10/3)*exp(- (2*kmrgd(1))/9 - 20/3)*(kmrgd(1)/500 + 3/50)*(kmrgd(3) - 1))/(243*(exp(- kmrgd(1)/9 - 10/3) - 1)^4) , 0 , exp(- kmrgd(1)/9 - 10/3)/(2250*(exp(- kmrgd(1)/9 - 10/3) - 1)^2) - (17*exp(- (10*kmrgd(1))/59 - 350/59))/34810 - (exp(- kmrgd(1)/9 - 10/3)*(kmrgd(1)/500 + 3/50))/(81*(exp(- kmrgd(1)/9 - 10/3) - 1)^2) + (2*exp(- (2*kmrgd(1))/9 - 20/3)*(kmrgd(1)/500 + 3/50))/(81*(exp(- kmrgd(1)/9 - 10/3) - 1)^3) , 0 ; (68*kmrgd(4)*exp(- (2*kmrgd(1))/21 - 24/7))/231525 - (11*exp(160/21 - (5*kmrgd(1))/63)*(kmrgd(4) - 1))/(294*(exp(160/21 - (5*kmrgd(1))/63) - 1)^2) + (11*exp(320/21 - (10*kmrgd(1))/63)*(kmrgd(4) - 1))/(147*(exp(160/21 - (5*kmrgd(1))/63) - 1)^3) + (125*exp(160/21 - (5*kmrgd(1))/63)*((99*kmrgd(1))/50 - 4752/25)*(kmrgd(4) - 1))/(250047*(exp(160/21 - (5*kmrgd(1))/63) - 1)^2) - (250*exp(320/21 - (10*kmrgd(1))/63)*((99*kmrgd(1))/50 - 4752/25)*(kmrgd(4) - 1))/(83349*(exp(160/21 - (5*kmrgd(1))/63) - 1)^3) + (250*exp(160/21 - (5*kmrgd(1))/63)*exp(320/21 - (10*kmrgd(1))/63)*((99*kmrgd(1))/50 - 4752/25)*(kmrgd(4) - 1))/(83349*(exp(160/21 - (5*kmrgd(1))/63) - 1)^4) , 0 , 0 , (11*exp(160/21 - (5*kmrgd(1))/63))/(35*(exp(160/21 - (5*kmrgd(1))/63) - 1)^2) - (34*exp(- (2*kmrgd(1))/21 - 24/7))/11025 - (25*exp(160/21 - (5*kmrgd(1))/63)*((99*kmrgd(1))/50 - 4752/25))/(3969*(exp(160/21 - (5*kmrgd(1))/63) - 1)^2) + (50*exp(320/21 - (10*kmrgd(1))/63)*((99*kmrgd(1))/50 - 4752/25))/(3969*(exp(160/21 - (5*kmrgd(1))/63) - 1)^3) ];
tens32=[ (250*exp(- (5*kmrgd(1))/57 - 175/114)*(kmrgd(1) - 58))/(1083*(exp(- (5*kmrgd(1))/57 - 175/114) + 1)^4) - (1000*exp(- (10*kmrgd(1))/57 - 175/57)*(kmrgd(1) - 58))/(1083*(exp(- (5*kmrgd(1))/57 - 175/114) + 1)^5) - (100*exp(- (5*kmrgd(1))/57 - 175/114))/(19*(exp(- (5*kmrgd(1))/57 - 175/114) + 1)^4) , 0 , 0 , 0 ; (47*exp(- (2*kmrgd(1))/11 - 62/11))/(1375*(exp(- (2*kmrgd(1))/11 - 62/11) - 1)^2) - exp(23/10 - kmrgd(1)/10)/40000 - (4*exp(- (2*kmrgd(1))/11 - 62/11)*((47*kmrgd(1))/500 + 1457/500))/(121*(exp(- (2*kmrgd(1))/11 - 62/11) - 1)^2) + (8*exp(- (4*kmrgd(1))/11 - 124/11)*((47*kmrgd(1))/500 + 1457/500))/(121*(exp(- (2*kmrgd(1))/11 - 62/11) - 1)^3) , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ];
tens33=[ 0 , 0 , -6*kmrgd(3)^2 , 0 ; 0 , 0 , 0 , 0 ; exp(- kmrgd(1)/9 - 10/3)/(2250*(exp(- kmrgd(1)/9 - 10/3) - 1)^2) - (17*exp(- (10*kmrgd(1))/59 - 350/59))/34810 - (exp(- kmrgd(1)/9 - 10/3)*(kmrgd(1)/500 + 3/50))/(81*(exp(- kmrgd(1)/9 - 10/3) - 1)^2) + (2*exp(- (2*kmrgd(1))/9 - 20/3)*(kmrgd(1)/500 + 3/50))/(81*(exp(- kmrgd(1)/9 - 10/3) - 1)^3) , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ];
tens34=[ 0 , 0 , 0 , -12 ; 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ; (11*exp(160/21 - (5*kmrgd(1))/63))/(35*(exp(160/21 - (5*kmrgd(1))/63) - 1)^2) - (34*exp(- (2*kmrgd(1))/21 - 24/7))/11025 - (25*exp(160/21 - (5*kmrgd(1))/63)*((99*kmrgd(1))/50 - 4752/25))/(3969*(exp(160/21 - (5*kmrgd(1))/63) - 1)^2) + (50*exp(320/21 - (10*kmrgd(1))/63)*((99*kmrgd(1))/50 - 4752/25))/(3969*(exp(160/21 - (5*kmrgd(1))/63) - 1)^3) , 0 , 0 , 0 ];
tens35=[ (250*exp(- (5*kmrgd(1))/57 - 175/114)*(kmrgd(1) - 58))/(1083*(exp(- (5*kmrgd(1))/57 - 175/114) + 1)^4) - (1000*exp(- (10*kmrgd(1))/57 - 175/57)*(kmrgd(1) - 58))/(1083*(exp(- (5*kmrgd(1))/57 - 175/114) + 1)^5) - (100*exp(- (5*kmrgd(1))/57 - 175/114))/(19*(exp(- (5*kmrgd(1))/57 - 175/114) + 1)^4) , 0 , 0 , 0 ; (47*exp(- (2*kmrgd(1))/11 - 62/11))/(1375*(exp(- (2*kmrgd(1))/11 - 62/11) - 1)^2) - exp(23/10 - kmrgd(1)/10)/40000 - (4*exp(- (2*kmrgd(1))/11 - 62/11)*((47*kmrgd(1))/500 + 1457/500))/(121*(exp(- (2*kmrgd(1))/11 - 62/11) - 1)^2) + (8*exp(- (4*kmrgd(1))/11 - 124/11)*((47*kmrgd(1))/500 + 1457/500))/(121*(exp(- (2*kmrgd(1))/11 - 62/11) - 1)^3) , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ];
tens36=[ 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ];
tens37=[ 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ];
tens38=[ 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ];
tens39=[ 0 , 0 , -6*kmrgd(3)^2 , 0 ; 0 , 0 , 0 , 0 ; exp(- kmrgd(1)/9 - 10/3)/(2250*(exp(- kmrgd(1)/9 - 10/3) - 1)^2) - (17*exp(- (10*kmrgd(1))/59 - 350/59))/34810 - (exp(- kmrgd(1)/9 - 10/3)*(kmrgd(1)/500 + 3/50))/(81*(exp(- kmrgd(1)/9 - 10/3) - 1)^2) + (2*exp(- (2*kmrgd(1))/9 - 20/3)*(kmrgd(1)/500 + 3/50))/(81*(exp(- kmrgd(1)/9 - 10/3) - 1)^3) , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ];
tens310=[ 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ];
tens311=[ -6*kmrgd(3)^2 , 0 , -12*kmrgd(3)*(kmrgd(1) + 80) , 0 ; 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ];
tens312=[ 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ];
tens313=[ 0 , 0 , 0 , -12 ; 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ; (11*exp(160/21 - (5*kmrgd(1))/63))/(35*(exp(160/21 - (5*kmrgd(1))/63) - 1)^2) - (34*exp(- (2*kmrgd(1))/21 - 24/7))/11025 - (25*exp(160/21 - (5*kmrgd(1))/63)*((99*kmrgd(1))/50 - 4752/25))/(3969*(exp(160/21 - (5*kmrgd(1))/63) - 1)^2) + (50*exp(320/21 - (10*kmrgd(1))/63)*((99*kmrgd(1))/50 - 4752/25))/(3969*(exp(160/21 - (5*kmrgd(1))/63) - 1)^3) , 0 , 0 , 0 ];
tens314=[ 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ];
tens315=[ 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ];
tens316=[ -12 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ; 0 , 0 , 0 , 0 ];
tens3(:,:,1,1) =tens31;
tens3(:,:,1,2) =tens32;
tens3(:,:,1,3) =tens33;
tens3(:,:,1,4) =tens34;
tens3(:,:,2,1) =tens35;
tens3(:,:,2,2) =tens36;
tens3(:,:,2,3) =tens37;
tens3(:,:,2,4) =tens38;
tens3(:,:,3,1) =tens39;
tens3(:,:,3,2) =tens310;
tens3(:,:,3,3) =tens311;
tens3(:,:,3,4) =tens312;
tens3(:,:,4,1) =tens313;
tens3(:,:,4,2) =tens314;
tens3(:,:,4,3) =tens315;
tens3(:,:,4,4) =tens316;
%---------------------------------------------------------------------------
function tens4  = der4(t,kmrgd,par_Cai,par_Iapp)
%---------------------------------------------------------------------------
function tens5  = der5(t,kmrgd,par_Cai,par_Iapp)