Did you read the help of polyfit ?
[p,S] = polyfit(x,y,n) returns the polynomial coefficients p and a structure S for use with polyval to obtain error estimates or predictions. Structure S contains fields R, df, and normr, for the triangular factor from a QR decomposition of the Vandermonde matrix of x, the degrees of freedom, and the norm of the residuals, respectively. If the data y are random, an estimate of the covariance matrix of p is (Rinv*Rinv')*normr^2/df, where Rinv is the inverse of R. If the errors in the data y are independent normal with constant variance, polyval produces error bounds that contain at least 50% of the predictions.
in your case, P must be the 2nd output argument.