f(0) = 1 exactly implies that p(1) = 1 exactly. So subtract that 1 from each y value.
y1 = y - 1;
Now, hand-waving time: subtract the p(1) from the polynomial, and factor the x. You get x * (p(6) * x^4 + p(5) * x^3 + p(4) * x^2 + p(3) * x + p(2)). So (wave, wave) divide your y1 by x,
y1divx = y1(2:end) ./ x(2:end); %remove the (0,1) pair
leaving (wave, wave) p(6) * x^4 + p(5) * x^3 + p(4) * x^2 + p(3) * x + p(2) to be fit by a 4th order polynomial:
p1divx = polyfit(x(2:end), y1div, 4);
y1divx_hat = polyval(p1divx, x);
y1_hat = y1divx_hat .* x + 1; %undo division by x and subtraction of 1
