Here is a straightforward way:
N = length(D);
DD = zeros(N);
for i=1:N
DD = DD + D^i;
end
DD = logical(DD);
DD(1:N+1:end) = false;
B = sum(DD,2);
Perhaps a link to explain the meaning of the powers of the adjacency matrix.
You can use the following code to visualize the resulting graph (note that the plot doesn't distinguish directed/undirected edges):
% circular layout
t = linspace(0,2*pi,N+1)'; t(end) = [];
xy = [cos(t) sin(t)];
% plot graph and label nodes
subplot(121), gplot(DD, xy, '-*')
text(xy(:,1), xy(:,2), num2str((1:N)'), 'BackgroundColor',[.4 .9 .5], ...
'VerticalAlign','bottom', 'HorizontalAlign','right')
axis square off
% adjacency matrix
subplot(122), spy(DD)
set(gca, 'XTick',1:N, 'YTick',1:N)
ylabel('from'), xlabel('to')