Design feature matrix of 8-neighborhood elements

I have a matrix (image) of size m x n that I am trying to rearrange differently. I want to design a feature matrix of size (m x n) x 9 where each row corresponds to a matrix element centered around its 8-neighborhood elements. My attempt is to loop sequentially through each element of the original matrix to extract the neighborhood values, however this method is too computationally heavy and takes too long to perform as the matrix size is very large. Is there anyway to do this cost-beneficially?

Attempt

M_feat = nan(size(img,1)*size(img,2), 9);
temp = zeros(size(img)+2);
temp(2:end-1,2:end-1) = double(img);

for i = 2:size(img,1)+1
       
    for j = 2:size(img,2)+1
        
        neighbors = temp(i-1:i+1, j-1:j+1); % read 3-by-3 mini-matrix
        neighbors = reshape(neighbors, 1, []); 
        M_feat((i-2)*size(img,1) + (j-1),:) = neighbors; % store row-wise
          
    end
    
end

Answer 1

[De Ar]

Got it!

new_img = zeros(size(img)+2);
new_img(2:end-1,2:end-1) = double(img);

% Image boundary coordinates without first/last row/column
inner_coord = [2 2; size(new_img,1)-1 size(new_img,2)-1];

% 9x2 matrix with 1 row for the relative shift to reach neighbors
[d2, d1] = meshgrid(-1:1, -1:1);
d = [d1(:) d2(:)];

% Cell array to store all 9 shifted images
temp = {};

for i = 1:size(d,1)
   % x-indices of the submatrix when shifted by d(i,1)
   coord_x = (inner_coord(1,1):inner_coord(2,1)) + d(i,1);
   % y-indices of the submatrix when shifted by d(i,2)
   coord_y = (inner_coord(1,2):inner_coord(2,2)) + d(i,2);
   % image matrices resulting from shift by d(i,)
   temp{i} = reshape(new_img(coord_x, coord_y), 1, []);
end

% Column-wise bind all 9 shifted images
M_feat = vertcat(temp{:}).';