Why my plot cannot show sigmoid graph? Whereas I already input to the function

Question

Hello I am confused with my code in MATLAB. I am use Runge-Kutta orde 4th. Why my plot cannot show sigmoid whereas I already input the sigmoid function into my code. Please somebody can explain how to fix it. This is my code.

clc;clear;
  
%input for time
t(1)=0;
dt=0.1; %time interval
t=0:dt:100; %time span

%input empty array
T=zeros(length(t),1); %empty array for t
M1=zeros(length(t),1); %empty array for M1
M2=zeros(length(t),1); %empty array for M2
M3=zeros(length(t),1); %empty array for M3
O=zeros(length(t),1); %empty array for O
P=zeros(length(t),1); %empty array for P

 %initial condition
    %M1 =10;
    %M2 = 0;
    %M3 = 0;
    %O =0;
    %P =0;

for j = 1:length((t))
    T(j+1)=T(j)+dt;
    M1(j+1)= M1(j)+1./(1+exp(-T(j))); % sigmoid eq
    
    k1M2 = dt*fRK4M2(M2(j),M1(j));
    k2M2 = dt*fRK4M2(M2(j)+k1M2/2,M1(j)+k1M2/2);
    k3M2 = dt*fRK4M2(M2(j)+k2M2/2,M1(j)+k2M2/2);
    k4M2 = dt*fRK4M2(M2(j)+k3M2,M1(j)+k3M2/2);
    M2(j+1) = M2(j)+1/6*(k1M2+2*k2M2+2*k3M2+k4M2);
    
    k1M3 = dt*fRK4M3(M1(j),M3(j));
    k2M3 = dt*fRK4M3(M3(j)+k1M3/2,M1(j)+k1M3/2);
    k3M3 = dt*fRK4M3(M3(j)+k2M3/2,M1(j)+k2M3/2);
    k4M3 = dt*fRK4M3(M3(j)+k3M3,M1(j)+k3M3);
    M3(j+1) = M3(j)+1/6*(k1M3+2*k2M3+2*k3M3+k4M3);
    
    k1O = dt*fRK4O(O(j),M1(j));
    k2O = dt*fRK4O(O(j)+k1O/2,M1(j)+k1O/2);
    k3O = dt*fRK4O(O(j)+k2O/2,M1(j)+k2O/2);
    k4O = dt*fRK4O(O(j)+k3O,M1(j)+k3O);
    O(j+1) = O(j)+1/6*(k1O+2*k2O+2*k3O+k4O);
    
    k1P = dt*fRK4P(P(j),M1(j));
    k2P = dt*fRK4P(P(j)+k1P/2,M1(j)+k1P/2);
    k3P = dt*fRK4P(P(j)+k2P/2,M1(j)+k2P/2);
    k4P = dt*fRK4P(P(j)+k3P,M1(j)+k3P/2);
    P(j+1) = P(j)+1/6*(k1P+2*k2P+2*k3P+k4P);
    
    k1M1= dt*fRK4M1(M1(j),M2(j),M3(j),O(j),P(j));
    k2M1= dt*fRK4M1(M2(j)+k1M2/2,M3(j)+k1M3/2,O(j)+k1O/2,P(j)+k1P/2,M1(j)+k1M1/2);
    k3M1= dt*fRK4M1(M2(j)+k2M2/2,M3(j)+k2M3/2,O(j)+k2O/2,P(j)+k2P/2,M1(j)+k2M1/2);
    k4M1= dt*fRK4M1(M2(j)+k3M2/2,M3(j)+k3M3/2,O(j)+k3O/2,P(j)+k3P/2,M1(j)+k3M1/2);
    M1(j+1) = M1(j)+(1/6*(k1M1+(2*k2M1)+(2*k3M1)+k4M1));
    end
  
figure;
plot (T,M1,'r','Linewidth',3)
xlabel('time')
ylabel('M1')

figure 
plot (T,M2,'b','Linewidth',3)
xlabel('time')
ylabel('M2')

figure 
plot (T,M3,'g','Linewidth',3)
xlabel('time')
ylabel('M3')

figure 
plot (T,O,'b','Linewidth',5)
xlabel('time')
ylabel('O')

%figure 
%plot (T,P,'r','Linewidth',5)
%xlabel('time')
%ylabel('P')

This is my function for M1, M2, M3, O, and P in separate file

**function M1 =fRK4M1(M1,M2,M3,O,P) ** 

  delta=50;
  K1= 10^-4;
  Ko=0.1;
  n=3;
  Oa=10;
  Pa=100;
  mu_1=10^-3;
  K2=5*10^-4;
  K3=10^-3;
  gamma=75;
 
  M1=(delta*M1*(1-(M1/gamma))-2*K1*M1*M1-M1*(K2.*M2)-((Oa-n)*K3*M1*M3)-((Pa-Oa)*Ko*M1*O)-(mu_1*M1)); 
    
**function M2 =fRK4M2(M2,M1)**

M1=10;
%M2=0;

mu_2=10^-3;
K1= 10^-4;
K2=5*10^-4;     

M2 = (K1*M1*M1)-(K2*M1*M2)-(mu_2*M2);

**function M3 =fRK4M3(M3,M1)**
%M1=10;
M2=0;
%M3=0;

mu_3=10^-3;
K2=5*10^-4;
K3=10^-3;

M3=(K2*M1*M2)-(K3*M1*M3)-(mu_3*M3);

**function O =fRK4O(O,M1)**

%M1=10;
M3=0;
%O=0;
Ko=0.1;
mu_o=10^-4;
K3=10^-3;

O = (K3*M1*M3)-(Ko*M1*O)-(mu_o*O);

**
function P =fRK4P(P,M1)**

%M1=10;
O=0;
%P=0;
Ko=0.1;
mu_p= 10^-5;

P  = (Ko*M1*O)-(mu_p*P);
 

my result like this my plot for M1,M2, and M3 and for O and P like this my plot for O and P