function [H,h] = Hilbert(n)
  H = zeros(n,n);
  h = zeros(n,1);
  for i=1:n
    for j=1:n
      H(i,j) = 1 / ( i+j-1);
      h(i) = h(i) + H(i,j);
    end
  end
