Faster calculation of a porwer sum

This function in my code is a bottleneck. Is it possible to make it run faster?

x and y are vectors of the same length and out is a matrix containing the power sum of x and y.

    function [out] = powerSum2(x,y,order)
        out = zeros(order+1);
        for ii = 0:order
            for jj = 0:order
                out(ii+1,jj+1) = sum(x.^ii.*y.^jj);
            end
        end
    end