ODE45 coupled variables

Hello,

I am trying to solve a coupled PDE using finite difference method. I have 2 variables, temperature (T) and concentration (Cbulk) that vary as a function of time (t) and distance (x). I have already solved this PDE for one variable (T) while keeping Cbulk constant.

The ODE call functions looks like

[t Tode] = ode15s('fd_rhs',tspan,T0); 

and theODE function looks like

% fd_rhs.m 
function ydot=fd_rhs(time,y) 
global nt R delH E A rho cp h k delr Cbulk s
for i=1:nt-1
    if i == 1
           dy(i) = k/(rho*cp*4.184*1000)*2*(s+1)*(y(i+1)-y(i))/(delr*delr)+delH**Cbulk;
    else 
dy(i) = k/(rho*cp*4.184*1000)*((y(i+1)-2*y(i) +y(i-1))/(delr*delr))+(s/R)*(y(i+1)-y(i-1))/(delr*delr)+delH**Cbulk; 
    end 
end
dy(nt)= (k/delr)/(h+(k/delr))*(k/(rho*cp*4.184*1000)*((y(nt)-2*y(nt-1) +y(nt-2))/(delr*delr)))+(s/R)*(y(nt)-y(nt-1))/(delr*delr)+delH**Cbulk;
ydot = dy';
end

My question is this: I want to incorporate a second variable (concentration) into my odefucntion fd_rhs. In the current solution using only one variable (temperature), the concentration is kept constant as Cbulk. Since fd_rhs has multiple ODE's for each grid point. how do I now incorporate a second variable into the same fd_rhs ? What should the function call and the fd_rhs look like ? Normally for standard coupled ODE's, there would be y(1) and y(2). But I already have y(1), y(2), y(3) ... y(nt) here for temperature at nt grid points.