Need help debugging so I can accurately portray graph of ODE system that is numerically solved

So I'm able to produce the following graph but I need help debugging so I can figure out why it doesn't look the way it should like in the following paper. I haver a feeling that it is due to the value of h, the step size, it may be too big. Everything else seems ok, but I need to be sure

Here is my graph:

Here is graph in research paper(only part a):

Here is my code:

function cells

%constants
epsilon=0.0001; %detachment rate
M=30; %total number of ligands
N=4; %max binding number
Ns1=20; %number of nuclei for first set
Ns2=10; %number of nuclei for second set

%computing time variable
T=500000; %length of integration
tstart=10^-4;
tend=10^5;
t=tstart;
n=500; %total time steps
logstart=log10(tstart);
logend=log10(tend);
lpm=[logstart:(logend-logstart)/n:logend]; %log plot mesh
iplot=1; %counter
tarray=zeros(1,1);
tarray(iplot)=t;
h=max(10^(-5),0.1*min(1,t));

%Numerical scheme for sigma<1/2
c=[Ns1;0;0;0;0]; %initial conditions
carray=zeros(N+1,1); %removed fixed number of columns for matrix in order to amend column vectors in loop dependent on t 
carray(:,iplot)=c; %adding initial condition to first slot of array 
m=M-sum([1 2 3 4 5]'.*carray(:,iplot));

while t<=tend %loop for top graph implemented with Classical RK4
    z=RK4(c,h,epsilon,m);
    iplot=iplot+1;
    carray=[carray z];
    m=M-sum([1 2 3 4 5]'.*carray(:,iplot));
    t=t+h;
    tarray=[tarray t];
    c=z;
    h=max(10^(-5),0.1*min(1,t));
end
length(tarray)
length(carray(2,:))
logt=log10(tarray);
figure(1)
plot(logt,carray(2,:),'g')
hold on
plot(logt,carray(3,:),'y')
plot(logt,carray(4,:),'b')
plot(logt,carray(5,:),'k')
hold off

%Classical RK4
function W4=RK4(c,h,epsilon,m) 
X1=c;
X2=c+h*(1/2)*rhs(X1,epsilon,m);
X3=c+h*(1/2)*rhs(X2,epsilon,m);
X4=c+h*rhs(X3,epsilon,m);
W4=c+h*((1/6)*rhs(X1,epsilon,m)+(1/3)*rhs(X2,epsilon,m)+(1/3)*rhs(X3,epsilon,m)+(1/6)*rhs(X4,epsilon,m));

%System of ODE's
function dcdt=rhs(c,epsilon,m)
dcdt=zeros(5,1);
dcdt(1)=-m*c(1)+epsilon*c(2);
dcdt(2)=-m*c(2)-epsilon*c(2)+m*c(1)+epsilon*c(3);
dcdt(3)=-m*c(3)-epsilon*c(3)+m*c(2)+epsilon*c(4);
dcdt(4)=-m*c(4)-epsilon*c(4)+m*c(3)+epsilon*c(5);
dcdt(5)=-epsilon*c(5)+m*c(4);

Answer 1

[Christopher Arreola]

I have fixed this, it was due to m