Can I use Runge-Kutta 4th order to solve the 3x3 co-dependent odes below?

• dNa/dt = -k0*exp(-ER/T)*Na
• dNb/dt = k0*exp(-ER/T)*Na
• dNw/dt = -k0*exp(-ER/T)*Na

Initial conditions: Na(0)=10 kmol, Nw(0)=20 kmol, Nb(0)=0 kmol

I have made some progress, but the RK4 doesn't seem to apply for Nb (?probably because the initial of Nb is 0?). My code is included below, along with the additional info needed, so please feel free to work on it, make corrections and/or suggestions for a different method.

Thank you in advance

clear
clc
close all

timestep=1.0;
t= 0:timestep:3600;
x_vec=zeros(3,length(t));
%%%arxikes sunthikes
xa0=10; %kmol
xb0=0; %kmol
xw0=20; %kmol
x_vec(:,1)=[xa0,xb0,xw0];
for idx= 1:length(t)-1
    idx/length(t)*3600;
    k1=Derivatives(x_vec(:,idx));
    k2=Derivatives(x_vec(:,idx)+k1*timestep/2);
    k3=Derivatives(x_vec(:,idx)+k2*timestep/2);
    k4=Derivatives(x_vec(:,idx)+k3*timestep);
    phi=(1/6)*(k1+2*k2+2*k3+k4);
    x_vec(:,idx+1)=x_vec(:,idx)+phi*timestep;
end
   x_vec(2,:)=x_vec(2,:)';
for a=1:3
figure()
plot(t,x_vec(a,:))
end

function dxdt_vec = Derivatives(x_vec)
k0=5.0E+10; %1/h
ER=9020;
T=298; %K

dxdt_vec=[-k0*exp(-ER/T):+k0*exp(-ER/T):-k0*exp(-ER/T)]*x_vec;