ODE property change at each step for cascade operation

Hello Matlab pros! I am working on a modelling thesis in Matlab and am experiencing a ton of difficulties trying to get the right numerical solution. My knowledge of Matlab is basic but I'll try and hopefully get the message through. Here goes: There is n coupled mixers in a row so as to analyze axial dispersion, defined by a system of differential equations for each. I want to change my number of mixers freely so as to best fit my experimental results, which should be possible for any number of mixers. Each system of differential equations gives a concentration and time dependance where the concentrations in the preceding mixer(s)simultaneously affect the solution in the next one and so on. The numerical solutions from the last mixer are then taken and integrated with CUMTRAPZ so as to provide a yield-time plot. I am using a loop for each mixer and so far there is a solution only for n=1. My problem is, how do I get the solutions from the previous loop for each time value back into the next ODE loop? Here is the code:

 n=input('number of mixers=');
  Y_d=0;        %the troublemaker-initial condition for n=1
  Z_0=[1 1 1.5049]; %initial conditions for ODE
  global theta_ex theta_ix Gx Kx rx   %fitting parameters
  theta_ex=0.0026;
  theta_ix=1.0833; 
  Kx=.8724;
  rx=.7;
  Gx=0.5049;
  Xt=0.74;
  taub=12.961;    %residence time
  for j=1:n
    [taux Z_1]=ode45(@rez,(0:0.01:taub),Z_0,[],Y_d,n,Xt);
    Y_d=Z_1(:,1); %THE "SHOULD BE" INPUT ARGUMENTS FOR NEXT LOOP
  end
  Y=Z_1(:,1)
  N_Y=cumtrapz(taux,Y);
  GR=Gx*rx/(1+Gx);
  F=GR*N_Y;
  e=1.48.*F;
  q=1.5485.*taux;

plot(q,e) 
function dZ_1dtaux = rez(taux,Z_1,Y_d,n,Xt)

global theta_ex theta_ix Gx Kx rx

if Z_1(2)>Xt
    Yrs=1;                       
else
    Yrs=Kx*Z_1(2);
end
if Z_1(2)~=Xt
    Jf=(Yrs-Z_1(1))/theta_ex;
elseif Z_1(2)==Xt && Z_1(1)<Kx*Xt
    Jf=(Yrs-Z_1(1))/theta_ex;
else
    Jf=0;
end
Js=(Z_1(3)-Z_1(2))/theta_ix;

dZ_1dtaux=[Jf-n*(Z_1(1)-Y_d)
           Js/rx-Gx*Jf
           -Js/(1-rx)];
end

Y_d should be values for Z_1(:,1) from the previous loop (mixer) and at the exact time step if possible. Has someone done something similiar? Any suggestions, formulations of the problem? Please help