UD1

n = size(a, 1);
S = {1};
D = Inf*(1:n);
P = cell(1,n);
while numel(S)
   p = S{end};
   S(end)=[];
   t = a(p(end),:);
   k = setdiff(find(t), p);
   for q = k
      r = [p q];
      S{end+1} = r;
      w = b(r(2:end),:)-b(r(1:end-1),:);
      d = sum(diag(w*w'));
      u=r(end);
      if d < D(u)
         D(u) = d;
         P{u} = r;
      end
   end
end

D(isinf(D))=0;
[i,j]=max(D);
c=P{j};