function xdot=glucfun4(t,x);
%  Glucose-Dependent Insulin Secretion full model
%  example call:  ode45('glucfun4',[0:80],[20;10;6;.3]);

global Vm Km Vs Ks Vmax2 K2 Vmax1 K1 L ko Io Go KI Ki tau n m

G=x(1,1);
I=x(2,1);
Gi=x(3,1);
J=x(4,1);
R1=(K1*Vmax1*(G-Gi)*I^n)/((K1+G)*(K1+Gi)*(Ki^n+I^n));
R2=(K2*Vmax2*(G*J^m-Gi))/((K2+G)*(K2+Gi));
Rm=(Vm*Gi)/(Km+Gi);
Rs=(Vs*(Rm^4+L^4))/(Ks^4+Rm^4+L^4);
Jinf=KI/(I+KI);
xdot(1,1)=-R1-R2-ko*(G-Go);
xdot(2,1)=Rs-ko*(I-Io);
xdot(3,1)=R1+R2-Rm;
xdot(4,1)=(Jinf-J)/tau;