function stoch2(mat,k,tlimit,m,groups);
%file with the matrices
%k=number of stages
%tlimit=time limit
%m - number of matrices per group
%groups-no- of groups of matrices
fid=fopen(mat,'r');
b=fscanf(fid,'%f',[k,k*m*groups]);
l=0;
for i=1:m*groups
    for j=1:k
       for p=1:k
           aa (j,p,i)= b(j,p+k*(i-1));
       end
   end
aa(:,:,i)=transpose(aa(:,:,i));   
end
fclose(fid);

for ia=1:groups %for each group of matrices
   for ja=1:m
      a(:,:,ja)=aa(:,:,(ja+(ia-1)*m));
   end       
   
%specify initial strucure vector
wvec=ones(k,1)/k;
w(:,1)=wvec;

%matrices=[a1(:) a2 (:)]
% generate sequence of environments
env=ceil(rand(1,tlimit)*m);

%initial condition
n=ones(k,1)/k;

%iteration
for t=1:tlimit
   c=a(:,:,env(t));
   %n=c*n;
   %r(t)=log(sum(n));
   %n=n/r(t);
%generate sequence of structure vectors
   wvec=c*wvec;
   r(t)=sum(wvec);
   wvec=wvec/r(t);
   w(:,t+1)=wvec;
end

loglambda=mean(r);

%specify initial reproductive value vector
vvec=ones(1,k)/k;
v(:,tlimit+1)=vvec';

% generate sequence of reproductive value vectors
for i=tlimit:-1:1 %i=1:tlimit 
   c=a(:,:,env(i));
   vvec=vvec*c;
   vvec=vvec/sum(vvec);
   v(:,i)=vvec';
end

% initialize the sensitivity and elasticity matrices
sensmat=zeros(k);
elasmat=zeros(k);

%calculate summations
for i=1:tlimit
   c=a(:,:,env(i));
   sensmat=sensmat+(v(:,i+1)*w(:,i)')/(r(i)*v(:,i+1)'*w(:,i+1));
   elasmat=elasmat+(v(:,i+1)*w(:,i)'.*c)/(r(i)*v(:,i+1)'*w(:,i+1));
end

% divide summations by timelimit
sensmat=sensmat/tlimit;
elasmat=elasmat/tlimit;
fid5=fopen('lambda','a');
fprintf(fid5,'%10.7f\t%s\n',loglambda,'mat');
fclose(fid5);
for i=1:k
   fid5=fopen('Stablestage.out','a');
   fprintf(fid5,'%10.7f\t',w(i,tlimit));
   if i==k
        fprintf(fid5,'\n');
   end  
	fclose(fid5);
   for j=1:k
       fid5=fopen('sen','a');
		 fprintf(fid5,'%10.7f\t%10.7f\t%10.7f\t%s\n',sensmat(i,j),i,j,'mat');
		 fclose(fid5);
		 fid5=fopen('ela','a');
		 fprintf(fid5,'%10.7f\t%10.7f\t%10.7f\t%s\n',elasmat(i,j),i,j,'mat');
       fclose(fid5);
    end
end
end