Psuedorandom variable causing problems during numerical solution and matrix concatanation

Monte Carlo Iterations of Numerical Solution to a System of Equations using a Psuedorandom variable

In the code below a psuedrandom distribution is applied to a single parameter (Y) which is explicitly considered in one equation of a dual equation system. The second equation implicitly considers the variable Y because it is linked to equation 1. Using a nested for loop a solution vector for both equations is constructed over 191 time steps. An outer loop is used to create a psuedorandom value for parameter Y and record this information in a 501-element row vector. This outer loop also builds a matrix containaing all 501 solution vectors for each equation. The resulting set of matrices should both be 501X191 however the matrix which explicitly considers the psuedorandom variable Y returns as a 1X1 array.

Why is this?

 %Initialize Time Set and Step Size

 h=0.01;        %h is the time step.
 tend=2;             %tend is end of approximation in days
 t=.1:h:tend;      %initialize time variable.

    % Construction of Matrix containing m solution vectors 

    % Generate Random Values of Y over a number of  
    % Monte Carlo iterations previously specified
 for j=1:(MCits-1),  

    Y(j+1)=(Y(j)*rand(1)/Y(j));   % Generate string of psuedorandom 
                                  % values of Y

    SUB(j,:)=vertcat(Substrate);  % construct matrix of approximate 
    MLSS(j,:)=vertcat(Biomass);   % solution vectors for all monte 
                                  % carlo iterations dim=(MCitsXt)

     clear Substrate;  %wipe out old variables.
     clear Biomass;

     Substrate(1)=S0;  %initial condition.
     Biomass(1)=X0;

    for i=1:length(t)-1,     %initialize for loop for time solution.

    %Calculate derivative (rate);
        Srate(i)=-((q*Substrate(i))/(K+Substrate(i)))*Biomass(i);  
        Xrate(i)=Y(j)*((q*Substrate(i))/(K+Substrate(i))-b-(k2h  
        /1.42))*Biomass(i);

     %Estimate new concentrations;
        Substrate(i+1)=Substrate(i)+h*Srate(i);  
        Biomass(i+1)=Biomass(i)+h*Xrate(i);
 end
 end

 S=size(SUB)  % answer returned 501 X 191

 *M=size(MLSS) % answer returned 1 X 1*

Note: The Xrate equation contains the psuedorandom variable Y. The Substrate equations depend on previous values of Xrate and Biomass but not explicitly on Y. The matrix is fully constructed for SUB (e.g. 501 rows of 191-element row vectors concatenated. An error is returned for Biomass which is explicitly dependent on Y.