Runge Kutta 2nd order 1st order ODE equations

Please help, those are the EDO´S i need to put in the graph, the graph should look like this:

I've been trying to fix it, but runge-kutta is a difficult method for me, thank you!

%RK de 2nd order 
%EDOS
% dx/dt = MU*xI - Kd*xI
%ds/dt = (-a/Yps)*(mu*xI);
%dp/dt = a*mu*xI;
h=25;
t=linspace(0,25,5e4);
%Initial indep. values
Kd = 0.0032;
a = 0.6212;
Yps = 0.5820;
Ki = 243.95
MU = 0.5557;
Ks = 0.0191;
p = 30.7600
xI =1.25;
sI =86.63;
pI =0;

%F_xt = @(x,t)  mu*xI - Kd*xI;  
%F_st = @(s,t) -a/Yps*mu*xI;  
%F_pt = @(p,t) a*mu*xI; 

for  i=1:length(t)-1
    mu(i) = MU*sI(i)/(Ks + sI(i)+ (sI(i)^2/Ki))*((1-pI(i)/p)^1);
    k1f1(i)= h*(mu(i)*xI(i) - Kd*xI(i));
    k2f1(i)= h*(mu(i)*(xI(i)+k1f1(i)) - Kd*(xI(i)+k1f1(i))); 
    xI(i+1) =xI(i) + 0.5*(k1f1(i) + k2f1(i));
    k1f2(i) = h*((-a/Yps)*mu(i)*xI(i));
    k2f2(i) = h*((-a/Yps)*mu(i)*(xI(i)));
    sI(i+1) =sI(i) + 0.5*(k1f2(i) + k2f2(i));
    k1f3(i) = h*(a*mu(i)*xI(i));
    k2f3(i) = h*(a*mu(i)*(xI(i)));
    pI(i+1) = pI(i) + 0.5*(k1f3(i) + k2f3(i));
   
end

 %grafica
 plot(t,xI,t,sI,t,pI);
legend('Cells','Subst','Product');

[Jan]

The code gets much easier, if you split the function to be integrated and the code to integrate. Do not run 3 different integrations, but combine the 3 components to a vector y =[x,s,p].

A problem is, that your ODE does not depend on the variables at all:

dx/dt = MU*xI - Kd*xI
ds/dt = (-a/Yps)*(mu*xI);
dp/dt = a*mu*xI;

These are all constants, which do not depend on x,s,p or t. What is mu in the the code?