Runge Kutta solving differential equations

I am trying to solve differential equations using runge kutta. I have two codes that I am working from . The first code I had an equation and dveloped runge kiutta from that equation. The second code I have four differential equations. I am unsure how to correctly incoroporate the k values from the first code into the second code to solve for k2, k3 and k4 (to which I have labled in the second code as e2,s2,p2,c2.....etc).

First code: %Define f(t,y)

f = @(t,y) (2 - exp(-4*t) - 2*y);

% Define Step Size and final time point (t_final)

h = 0.001;

t_final = 5;

t = 0:h:t_final;

y = zeros(1,numel(t));

y(1) = 1; % y0

% You know the value a t = 0, thats why you'll start: t = h i.e. i = 2

for i = 2:numel(t)

    k1 = h*f(t(i-1),y(i-1));

    k2 = h*f(t(i-1)+h/2, y(i-1)+k1/2);

    k3 = h*f(t(i-1)+h/2, y(i-1)+k2/2);

    k4 = h*f(t(i-1)+h, y(i-1)+k3);

     y(i) = y(i-1) + (k1+2*k2+2*k3+k4)/6;

     disp([t(i) y(i)]);
end

plot(t,y);

Second code: %define your f(t,y)

h = 0.1; % Define Step Size

t = 0:h:200; %define your range of time points (range of t values)

S(1)=100;

E(1)=10;

C(1)=0;

P(1)=0;

k1=0.005;

k2=k1;

k3=0.1;

% You know the value a t = 0, thats why you'll state with t = h i.e. i = 2 %Write for loop for i = 1:numel(t)-1

   e1=(k2*C(i)+k3*C(i)-k1*E(i)*S(i)); 

   s1=(k2*C(i)-k1*E(i)*S(i));

   c1= k1*S(i) * E(i) -k2*C(i) -k3*C(i) -k2*C(i);

   p1=(k3*C(i));

    e2=

    s2=

    c2=

    p2=

    e3=

    s3=

    c3=

    p3=

    e4=

    s4=

    c4=

     p4=

   E(i+1) = E(i) + h*e1;

    S(i+1) = S(i) + h*s1;

    C(i+1) = C(i) + h*c1;

    P(i+1)= P(i)+ h*p1;

end % plot each graph plot(t,E)

hold on

plot(t,S)

plot(t,C)

plot(t,P)