function Y = mds(D)
n = size(D,1);
A = -0.5*D;
H = eye(n)-(1/n)*ones(n);
B = H*A*H;
[V,L] = eig(B); % V: eigenvector. L:eigenvalue
L = real(diag(L)); % stored as real number
[L,I] = sort(L,'descend'); % rearrange in decreasing order, record the original position by I
Y = zeros(n);
for i = 1:3
Y(:,i) = V(:,I(i))*(L(i)^0.5);
end
hold on;
for i = 1:n
plot3(Y(i,1),Y(i,2),Y(i,3),'.', 'MarkerSize', 16, 'Color',[0 1-i/n i/n]);
end