Finding appropriate step size for set of ODE's numerically solved without an ODE solver

So I'm trying to implement the following code to replicate this graph(parts a. and c. only)

Here is my code, the time step is increasing as we move from from 10^-4 to 10^5. So the step size, h, should be increasing in size as it goes along. I would think I would need an adaptive time step, but don't know how to implement unless I use the fancier Runge-Kutta Fehlberg 4/5 pair. I did post earlier but was left hanging with regards to this question, I'm trying to do this without using an ODE solver like ODE15 or ODE45

function cells

%constants
epsilon=0.0001;
M=30;
N=4;%double check if it plays role in code, mentioned in legend for graph
Ns1=20;
Ns2=10;
m=M;%constraint condition

%computing time variable
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,n+1);
tarray(iplot)=t;
%h=(10^-4);

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

while t<tend %loop for top graph implemented with Classical RK4
    z=RK4(c,h,epsilon,m);
    tm=10^(lpm(iplot));
    if t>tm
        iplot=iplot+1;
        carray(:,iplot)=z;
        tarray(iplot)=t;
    end
    c=z;
    t=t+h;
end

logt=log10(tarray);
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
c=[Ns2,0,0,0,0];
carray=zeros(5,n+1);
carray(:,1)=c;

while t<tend%loop for bottom graph implemented with Classical RK4
    z=RK4(c,h,epsilon,m);
    tm=10^(lpm(iplot));
    if t>tm
        iplot=iplot+1;
        carray(:,iplot)=z;
        tarray(iplot)=t;
    end
    c=z;
    t=t+h;
end

logt=log10(tarray);
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);