function J = van(p,x)
% VAN Vandermonde matrix
%     J = van(p,x) computes p+1 rows
%     of the Vandermonde of the vector x.
n1 = max(size(x));
J = ones(n1,1);
for j = 1:p
   J = [x.^j J];
end

