Reproduce this figure using matlab for system of ODE's solved numerically

Question

So I believe that I am working with the Becker-Doring equations or the spatially uniform mass-action equations, either way, this was taken from a research paper and I'm trying to reproduce the results using a numerical method of my choice to solve the system of ODE's and then reproduce the figures. I will most likely use Runge-Kutta Method of order 4. Just hoping someone on here has the skills and knows exactly how to reproduce this figure in matlab cause my skills are beginner level.

Here is my code so far, basic

function cells
epsilon=0.0001;
M=30;
N=4;
Ns1=20;
Ns2=10;
X=[.0001 .001 .01 .1 1 10 100 1000 10000 100000];
Y=log10(X);
c=zeros(4,1);
m=M;
for t=1:4
    m=m-t*c(t,1);
    t=t+1;
end
for t=1:4
    z=c(t,:)+h*euler(c(t,:),epsilon,m); 
    c=[c z];
    t=t+1;
end
figure(1)
plot(Y,c(1,:),'g');
hold on
plot(Y,c(2,:),'r');
plot(Y,c(3,:),'b');
plot(Y,c(4,:),'k')
hold off
function y=euler(z,epsilon,m)
y(1,1)=-m*z(1,:)-epsilon*z(1,:)+epsilon*z(2,:);
y(2,1)=-m*z(2,:)-epsilon*z(2,:)+m*z(1,:)+epsilon*z(3,:);
y(3,1)=-m*z(3,:)-epsilon*z(3,:)++m*z(2,:)+epsilon*z(4,:);
y(4,1)=-m*z(4,:)-epsilon*z(4,:)+m*z(3,:);

Here is the link to the entire paper

Update on code: Overall, I think I'm getting closer to figuring this out. I switched from Euler to Runge-Kutta for my numerical scheme.

function cells

%constants
epsilon=0.0001;
M=30;
N=4;
Ns1=20;
Ns2=10;
m=M;%constraint condition

%computing time variable
tstart=10^-4;%Need to implement adaptive time step
tend=10^5;
t=tstart;
h=10^-4;%uniform step size
n=round((tend-tstart)/h);
tspan=[10^-4 10^5];%replace after time step is figured out
logt=log10(tspan);

%Numerical scheme for sigma<1/2
c0=[Ns1,0,0,0,0];%initial conditions
carray=zeros(5,n+1);
carray(:,1)=c0;%adding initial condition to array

for i=2:length(logt)%loop for top graph implemented with Classical RK4
    z=RK4(c,h,epsilon,m); 
    carray(:,i)=z;
    c=z;
    i=i+1;
    t=t+h;
end

figure(1)
tiledlayout(2,1) %initiating display of multiple axes in one figure
nexttile%top plot
plot(logt,carray(2,:),'g')
hold on
plot(logt,carray(3,:),'y')
plot(logt,carray(4,:),'b')
plot(logt,carray(5,:),'k')
hold off

%%Numerical scheme for sigma>1/2
c0=[Ns2,0,0,0,0];
carray=zeros(5,n+1);
carray(:,1)=c0;

for i=2:length(logt)%loop for bottom graph implemented with Classical RK4
    z=RK4(c,h,epsilon,m); 
    carray(:,i)=z;
    c=z;
    i=i+1;
    t=t+h;
end

nexttile%bottom plot
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(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);