function repressilator

   global parameters % Share with subroutines/functions

   %parameters = [0,50,0.2,2]; % [alpha0,alpha,beta,n] for Fig. 4.2
   %parameters = [0,1,0.2,2]; % [alpha0,alpha,beta,n] for Fig. 4.3
   parameters = [1,50,0.2,2]; % [alpha0,alpha,beta,n] for Fig. 4.3
   
   u0 = [0.2,0.1,0.3,0.1,0.4,0.5]; % Initial conditions

   int_time = [0 300]; % Integration time.
   options = odeset('AbsTol',1e-8,'RelTol',1e-8); % Control accuracy of solver
   [t,u] = ode45(@f,int_time,u0,options);

   figure(1) 
   subplot(3,1,1)
      plot(t,u(:,1),'b',t,u(:,2),'k');
      xlabel('time')
      ylabel('m1 blue, p1 black')
   subplot(3,1,2)
      plot(t,u(:,3),'b',t,u(:,4),'k');
      xlabel('time')
      ylabel('m2 blue, p2 black')
   subplot(3,1,3)  
      plot(t,u(:,5),'b',t,u(:,6),'k');
      xlabel('time')
      ylabel('m3 blue, p3 black')
   figure(2)      
      plot(t,u(:,1),'b',t,u(:,3),'g',t,u(:,5),'k');
      xlabel('time')
      ylabel('m1 blue, m2 green, m3 black');
end

function dudt = f(t,u)

   global parameters
   dudt = zeros(6,1); % Create a column vector.

   a0 = parameters(1);
   a  = parameters(2);
   b  = parameters(3);
   n  = parameters(4);

   m1=u(1); p1=u(2);
   m2=u(3); p2=u(4);
   m3=u(5); p3=u(6);

   dudt(1) = -m1 + a/(1+p3^n) + a0;
   dudt(2) = -b*(p1-m1);
   dudt(3) = -m2 + a/(1+p1^n) + a0;
   dudt(4) = -b*(p2-m2);
   dudt(5) = -m3 + a/(1+p2^n) + a0;
   dudt(6) = -b*(p3-m3);	

end
