I don't understand the rules of evolution of the board that you described. However, an efficient way to gather all the 8-neighbor groups for subsequent analysis is as follows. Here, I've used a 4x4 board size for illustration,
n=4; %board size
[I,J,dI,dJ]=ndgrid(1:n,1:n,-1:1,-1:1);
I=reshape(I+dI,[],9); clear dI
J=reshape(J+dJ,[],9); clear dJ
valid=(I>=1 & I<=n) & (J>=1 & J<=n);
I(~valid)=nan;
J(~valid)=nan;
lookup=sub2ind([n,n],I,J);
lookup(~valid)=1;
lookup(:,5)=[];
valid(:,5)=[];
mask=double(valid); mask(~valid)=nan;
The key results of this are lookup and mask. You would compute them only once and pass them to your function for repeated use. You use them as follows, to form a matrix whose rows are the neighbors of successive locations on the board. The NaNs indicate nieghbors outside the boundaries of the board.
board=rand(n),
neighbors = board(lookup).*mask
