%this routine solves a system of odes

clear
global vm vp km k3 k1 k2 J K A

%specify parameter values
vm = 1;
vp = 0.5;
km = .1;
k3 = .1;
k1 = 10;
k2 = 0.03;
J = 0.05;
K = 200;
A = 0.1;

P = [ 0:.05:4];
q = 2./(1+sqrt(1+8*K*P));
P1 = q.*P;
P2 = ((1-q)/2).*P;

M1  = vm./((1+(P2/A).^2) * km);

M2 =((k1*P1+2*k2*P2)./(J+P)+k3*P)/vp;

%specify the output points
tspan = [0:.1:100];

%specify initial data

y0 = [0,0];



[T,Y] = ode23('circadian_deRHS',tspan, y0);  
figure(1)
plot(T,Y(:,1),'r',T,Y(:,2),'b','linewidth',2)
xlabel('time (hours)','fontsize',16)
legend('mRNA','Per')


figure(2)
plot(M1,P,'b--',M2,P,'r--',Y(:,1),Y(:,2),'linewidth',2)
xlabel('mRNA','fontsize',16)
ylabel('Per','fontsize',16)
text(5,.4,'dM/dt = 0','fontsize',16)
text(1.5,3.5,'dP/dt = 0','fontsize',16)
