I understand I can get all linear trends at once using the backslash operator. But that doesn't give p-values.
First, it might not be a good idea to use backslash straightforwardly. Polynomial fitting has certain numerically sensitive properties. POLYFIT uses a QR decomposition:
[Q,R] = qr(V,0);
ws = warning('off','all');
p = R\(Q'*y);
if size(R,2) > size(R,1)
if nargout > 2
But anyway, once you've found a solution to the system V*x=Y, it seems to me that you could vectorize the p-values by analyzing the residuals,