global a b c d e f;
a=0.2;b=1;c=0.5;d=0.1;e=0.3;f=15;
%initial conditions
S0=15;X0=0.05;P0=0;tend=100;
y0=[X0 P0 S0];

[t,y]=ode45('myfun',[0 tend],y0);
%results
X=y(:1);P=y(:2);S=(:3);

plot(t,S);