clear;
close;
To=8e-6;
Yo=20e-6;
Bo=1.7e-6;

options = odeset('RelTol',1e-9,'AbsTol',[1e-9 1e-9 1e-9 1e-9 1e-9]);

[t y]=ode23('spirofunc',[0 80],[4e-6 4e-6 0e-6 0.5e-6 10e-6],options);

Ptot=To-y(:,1)-y(:,2);
Bp=y(:,4);
Yp=y(:,5);
metlevel=1-(y(:,1)+y(:,3))/To;
phoslevel=1-(y(:,1)+y(:,2))/To;
subplot(2,2,1)
plot(t,phoslevel,'bx');
axis([10 80 0 0.1]);
title('Phosphorylation Level');
subplot(2,2,2)
plot(t,metlevel,'rx');
title('Methylation Level');
axis([10 80 0 1]);
subplot(2,2,3)
plot(t,Bp/Bo,'gx');
axis([10 80 0 1]);
title('Bp/Btot');
subplot(2,2,4)
plot(t,Yp/Yo,'yx');
axis([10 80 0 1]);
title('Yp/Ytot');