Assuming A is the transposed of the adjacent matrix (n x n).
j in (1:n) is the source node number.
This method should be faster than computing A^k proposed by Matt
x = zeros(n,1);
x(j) = 1;
x = A*x;
i = find(x)';
fprintf('%d-neighbor of %d is %s\n', k, j, mat2str(i))
If sparse form of the adjadcent matrix is readily available, it could be preferable to use it.