How do you add demographic stochasticity to a stage structure population model?

Hi, I have a 3x3 transition matrix and need to incorporate demographic stochasticity. I tried using binornd(n,p)/n to get a random proportion of stage 1 individuals transitioning each year, but it doesn't appear to be working. It doesn't look like the plot is calling this every year.

Thanks in adance

if true
  % code
end
% The demographic parameters
m33 = .98 ;
P21_g = 0.908 ;  P21_b = 0.908  ; P21 = 0.908 ;
P22_g = 0.944*0.75 ; P22_b = 0.944*0.75 ; P22 = 0.944*0.75 ;
P32_g = 0.944*0.25 ; P32_b = 0.944*0.25 ; P32 = 0.944*0.25 ;
P33_g = 0.97 ; P33_b = 0.97*m33 ; P33 = 0.97 ;
F3 = 0.5*P33*0.044 ; % 0.5 probability of calf being female * survivability of mature * survivability of mother * fecundity (0.044, from pace)
n1 = 16 ; n2 = 30 ; n3 = 140 ; % the initial condition
N0 = [n1 n2 n3]' ; % note that N0 is written as a column vector to allow matrix calculations
% The three projection matrices i am considering for environmental
% stochasticity (A_g = good year, no mortality event; A_b = bad year,
% mortality event)
A_g = [0                      0      F3      ;
    binornd(16, P21_g)/16  P22_g   0       ;
     0                     P32_g  P33_g]   ;    
A_b = [0                      0      F3      ;
    binornd(16, P21_b)/16  P22_b   0       ;
     0                     P32_b  P33_b]   ;
A = [0                    0    F3      ;
  binornd(16, P21)/16  P22   0       ;
   0                   P32  P33]     ;
% ===== Plotting the population time series and population growth rates for 40 years 
% ===== assuming a stable (deterministic) and average environment
N_ts             = zeros(3, 21) ; % a matrix that will contain the number of individuals in each age class (columns) over time (rows)
N_ts_total       = zeros(1, 21) ; % a vector that will contain the total population size as a function of timee 
lambda           = zeros(1, 21) ; % a vector that will contain the annual population growth rates as a function of time
N_ts(:, 1)       = N0           ; % the initial condition is written into the first column of N_ts 
N_ts_total(1, 1) = sum(N0)      ; % the initial total population size is the sum of the initial abundances in all age classes 
% Project the population 40 years into the future using matrix projections
for i = 2 : 21
  N_ts(:, i)       = A * N_ts(:, i-1) ; % update the population vector by one year
  N_ts_total(1, i) = sum(N_ts(:, i))  ; % sum all individuals in the four age classes in a year to get total pop size 
  lambda(1, i)     = N_ts_total(1, i) / N_ts_total(1, i-1) ; % annual pop growth rate is the the population size in one year, divided by the size in previous year
end
% time series - stage class
figure(1)
plot(0 : 20, N_ts)
hold on
legend ('calves','imm. females','mat. females')
xlabel('\it_t')    
ylabel('abundance')