\(a9 = 20; am = 4; lead = 1.05; 
p[0] = 1 - 0.687*lead; 
    p[1] = 1 - 0.248*lead; p[2] = 1 - 0.215*lead; p[3] = 1 - 0.066*lead; 

    Do[p[a] = p[3], {a, 4, a9}]; 
Do[m[a] = 0. , {a, \(-20\), am - 1}]; 

    Do[m[a] = m1, {a, am, a9}]; 
Do[l[i, j] = 0, {i, 1, a9}, {j, 1, a9}]; 

    \(Do[l[1, j] = m[j]\ p[0], {j, 1, a9}]; \)\)\(Do[l[j + 1, j] = p[j], {j, 1, a9}]; 
\(l[a9, a9] = p[j]; 
    \) (*NOT\ die\ at\ 20\ yrs\ old*) ; 

    ll = Table[l[i, j], {i, 1, a9}, {j, 1, a9}]; 

    m10 = \((0.268 + 2*0.261)\)/2; m11 = m10/4; 

    Print[{m10, 
        \(Eigenvalues[ll /. m1 -> m10]\)[[1
          ]], 
        \(Eigenvalues[ll /. m1 -> m11]\)[[1
          ]]}]; \)

n0=Eigenvectors[ll/.m1->m10][[1]];
ns[0]=Sum[n0[[a]],{a,1,a9}];
n0=n0/Sum[n0[[a]],{a,2,a9}]*1500;
ns[0]=1500;

MatrixForm[N[Table[l[i,j],{i,1,6},{j,1,6}]/.m1->m10,3]]
N[n0,2]

k9=100;Do[n[0]=n0;f1=0.75;t9=100;
Do[n[t+1]=N[ll/.m1->If[Random[]<f1,m10,m11]].n[t];
ns[t+1]=Sum[n[t+1][[a]],{a,2,a9}],{t,0,t9-1}];
gr[k]=ListPlot[Table[ns[t],{t,0,t9}],PlotJoined->True,DisplayFunction->Identity];n9[k]=ns[t9],{k,1,k9}];
Show[Table[gr[k],{k,1,k9}],PlotRange->{{0,100},{0,Automatic}},DisplayFunction->$DisplayFunction];

a9=100;
i9=200;(*initial population size*);
f1=0.75;(*freq of good year*);t9=100;(*end of years*);
k9=1000;(*# of simulations*);
Do[risk[t]=0;nm[t]=If[t==0,i9* k9,0],{t,0,t9}];(*risk=extinction risk until yr t, nm=mean population size in yr t * # runs*);
Do[Do[h[i]=0,{i,1,300}];(*state of individual i*);
Do[h[i]=Mod[i-1,3]+2,{i,1,i9}];(*the absolute quantity in the left hand means age, negative number means male*);
Do[h[i]=-h[i],{i,1,i9/2}];
ns[0]=i9;(*# of individuals in yr 0*);
Do[ns9=ns[t-1];
Do[If[Random[]>p[Min[4,Abs[h[i]]]]||Abs[i]==a9,{h[i]=h[ns9];h[ns9]=0;ns9=ns9-1}],{i,1,ns9}];(*death prob.*);
female=If[ns9>0,Sum[If[h[i]>3,1,0],{i,1,ns9}],0];
(*# of male adults*);
male=If[ns9>0,Sum[If[h[i]<-3,1,0],{i,1,ns9}],0];
cpl=Min[male,female];(**);
m0=If[Random[]<f1,m10*2,m11*2];cl=0;cp=0;(*# of fledgings per mother, # of total fledgings*);
Do[If[cp<cpl&&Random[]<m0&&h[i]>3,{cl=cl+2,cp=cp+1}],{i,1,ns9}];cl1=cl;(*make pair*);
Do[If[Random[]>p[0],cl=cl-1],{i,1,cl1}];(*death until 1 year old*);
Do[h[i]=Sign[h[i]]*(Abs[h[i]]+1),{i,1,ns9}];
Do[h[ns9+i]=If[Random[]<0.5,1,-1],{i,1,cl}];(*sex ratio*);
ns[t]=ns9+cl;If[ns[t]==0,risk[t]=risk[t]+1];nm[t]=nm[t]+ns[t],{t,1,t9}];
gr[k]=ListPlot[Table[{t,ns[t]},{t,0,t9}],PlotJoined->True,DisplayFunction->Identity],{k,1,k9}]

hh=Table[h[i],{i,1,ns[t9]}];
{ns[t9],Sort[hh]}

gr0=Show[Table[gr[k],{k,1,100}],PlotRange->{{0,t9},{0,300}},DisplayFunction->$DisplayFunction];

gr1=ListPlot[Table[risk[t]/k9,{t,0,t9}],PlotJoined->True];
gr2=ListPlot[Table[nm[t]/k9,{t,0,t9}],PlotJoined->True];