syms p;     %% prime in question
syms times; %% (number of iterations)-1
p=5;
times=3;
pack;



mat = repmat(0,[p^times p^times (times+1]));
mat(1,1,1)=1;

for reps=1:times
   for row=1:p^(reps-1)
       for col=1:p^(reps-1)
           
           %% reproduce into the next matrix
           %% starting with the (R,C)th entry and going p places 
           % in each direction in the new matrix
           
           % fill in the LEFT COLUMN of tile
           C=(col-1)*p+1;
           for R=(row-1)*p+1:(row-1)*p+p
                   
               lambda = mod((1/p)*(NCHOOSEK(p+mod(R-1,p),p)-1), p);
               val = mat(row,col,reps) * (1+p*lambda*mod(row,p));
               mat(R,C,reps+1) = mod(val,p^2);
               pack;
               
           end
           
           % fill in the TOP ROW of tile
           R=(col-1)*p+1;
           for C=(col-1)*p+1:(col-1)*p+p
               
              lambda = mod((1/p)*(NCHOOSEK(p+mod(C-1,p),p)-1), p);
               val = mat(row,col,reps) * (1+p*lambda*mod(row,p));
               mat(R,C,reps+1) = mod(val,p^2);
               pack;
               
           end        
              
           %% fil in rest of tile according to assembly rule
           for R=(row-1)*p+2:(row-1)*p+p 
             for C=(col-1)*p+2:(col-1)*p+p
          
                 mat(R,C,reps+1) = mod(mat(R,C-1,reps+1)+mat(R-1,C,reps+1), p^2);
                 pack;
                 
             end 
         end
                   

           
                  
           
     end
 end
end
