A = [x(:).^2, x(:), ones(numel(x),1)];
b = y;
coefficients = (A\b).';
This code is fine with x being either a row vector or column vector. The code expects that you are doing the quadratic fitting over each column of y independently (not over the rows.)
The array of coefficients that results will be N x 3 where N is the number of columns of y. The order of the coefficients will be the x^2 coefficient, then the x coefficient, then the constant.