function metabol();
x0=[0.20000000D+01    0.30000000D+02    0.12482071D+01    0.49375896D+02    0.49375896D+02    0.98638884D+02    0.13611157D+01
    0.30000000D+02    0.25000000D+01    0.12482071D+01    0.49375896D+02    0.49375896D+02    0.98638884D+02    0.13611157D+01
    0.20000000D+02    0.30000000D+02    0.12482071D+01    0.49375896D+02    0.49375896D+02    0.98638884D+02    0.13611157D+01
    0.10000000D+00    0.30000000D+02    0.12482071D+01    0.49375896D+02    0.49375896D+02    0.98638884D+02    0.13611157D+01
    0.10000000D+00    0.30000000D+02    0.13240875D+01    0.17474971D+00    0.98501163D+02    0.98579949D+02    0.14200506D+01
    0.10000000D+01    0.20000000D+01    0.14762975D+02    0.14762975D+02    0.70474050D+02    0.75579373D+01    0.92442063D+02
    0.10000000D+01    0.10000000D+01    0.98187001D+02    0.90649930D+00    0.90649930D+00    0.53093874D+00    0.99469061D+02
    0.20000000D+01    0.20000000D+01    0.42660956D+02    0.42660956D+02    0.14678088D+02    0.13937337D+01    0.98606266D+02
    0.10000000D+02    0.50000000D+00    0.98189998D+02    0.90650088D+00    0.90650088D+00    0.53092125D+00    0.99469979D+02
    0.50000000D+00    0.10000000D+02    0.98187001D+02    0.90649930D+00    0.90649930D+00    0.53093874D+00    0.99469061D+02
    0.15000000D+01    0.30000000D+02    0.12482071D+01    0.49375896D+02    0.49375896D+02    0.98638884D+02    0.13611157D+01
    0.50000000D+00    0.25000000D+01    0.13237984D+01    0.98150407D+02    0.52579430D+00    0.98580176D+02    0.14198240D+01];
tspan = [0,25,50,75,100,125,150];
for iexp=1:12
    [t,p] = ode15s(@ode_metabolic, tspan, x0(iexp,:));
    figure(iexp);
    plot(t,p(:,3),'ro-');
    hold on;
    plot(t,p(:,4),'bs-');
    plot(t,p(:,5),'gd-');
    plot(t,p(:,6),'k^-');
    plot(t,p(:,7),'m+-');
    legend('X_3','X_4','X_5','X_6','X_7');
end


function dx = ode_metabolic(t,X);
dx = zeros(7,1);
v1=0.10E+01*X( 4)/( 0.50E+01+X( 4));
v2=0.10E+01*X( 5)/( 0.50E+01+X( 5));    
v3=0.50E+01*X( 3)/((X( 3)+ 0.50E+01)*(1+X( 1)/ 0.10E+01));
v4=0.50E+01*X( 3)/((X( 3)+ 0.50E+01)*(1+X( 2)/ 0.10E+01));   
v5=0.10E+01*X( 6)/( 0.50E+01+X( 6));  
v6=0.10E+02*X( 7)/((X( 7)+ 0.50E+01)*(1+X( 3)/ 0.10E+01));  
dx(1)=0; % independent
dx(2)=0; % independent
dx(3)=v1+v2-v3-v4;
dx(4)=-v1+v3;
dx(5)=-v2+v4;
dx(6)=-v5+v6;
dx(7)=v5-v6;

