function ss_15genes();
n=15;
alpha=ones(n,1);
beta=ones(n,1);
g=zeros(n,n);
h=zeros(n,n);
g(1,14)=-0.1;
g(3,12)=-0.2;
g(5,1)=1.0;
g(6,1)=1.0;
g(7,2)=0.5;
g(7,3)=0.4;
g(8,4)=0.2;
g(9,5)=1.0 ;
g(9,6)=-0.1;
g(10,7)=0.3;
g(11,4)=0.4;
g(11,7)=-0.2;
g(12,13)=0.5;
g(13,8)=0.6;
g(14,9)=1.0;
g(14,15)=-0.2;
g(15,10)=0.2;
for i=1:n
    h(i,i)=1;
end
p0=[0.19002585702943E+01    0.11578260956854E+00    0.16762368901048E+01 ...
    0.99310489940621E+00    0.15896421604019E+01    0.11655833631225E+01 ...
    0.83074972088645E+00    0.13666646486770E+01    0.12170807224480E+01 ...
    0.24209422607283E+00    0.87958171389276E+00    0.27401978341524E+00 ...
    0.59481135791852E+00    0.38778635996934E+00    0.14330847945001E+01];
tspan = [0:0.2:2];
[t,p] = ode15s(@ode_ssys, tspan, p0,[],alpha,beta,g,h);
figure(1);
plot(t,p(:,1),'rs-');; 
hold on;
plot(t,p(:,2),'gd-');
plot(t,p(:,3),'ko-');
plot(t,p(:,4),'mx-');
plot(t,p(:,5),'b*-');
legend('X1','X2','X3','X4','X5');
figure(2);
plot(t,p(:,6),'rs-');; 
hold on;
plot(t,p(:,7),'gd-');
plot(t,p(:,8),'ko-');
plot(t,p(:,9),'mx-');
plot(t,p(:,10),'b*-');
legend('X6','X7','X8','X9','X10');
figure(3);
plot(t,p(:,11),'rs-');; 
hold on;
plot(t,p(:,12),'gd-');
plot(t,p(:,13),'ko-');
plot(t,p(:,14),'mx-');
plot(t,p(:,15),'b*-');
legend('X11','X12','X13','X14','X15');
xlabel('Time')
xlabel('Concentration')
title('A genetic network');


function dp = ode_ssys(t,X,alpha,beta,g,h);
dp = zeros(15,1);   
for i=1:15   
  vplus=alpha(i); 
  for j=1:15
     vplus=vplus*X(j)^g(i,j);
  end
  vminus=beta(i);  
  for j=1:15
     vminus=vminus*X(j)^h(i,j);
  end
  dp(i)=vplus-vminus;
end
