A power model doesn't fit these data worth a hooey at the lower end...try
plot(f,x,y,'*')
set(gca,'yscale','log')
ylim([y(1)/2 y(end)*2])
to see the problem; the curvature is positive for the fitted curve but negative in that region for the actual data. It appears either a 'power1' or linear in both log(x) log(y) would be a better fit of these data than 'power2'.
Addendum
b=polyfit(log(x),log(y),1); % fit log-log
yh=exp(polyval(b,log([x(1) x(end)]))); % Evaluate at end points
hold on
plot([x(1) x(end)],yh,'c') % add to the previous plot
You'll see this is quite a lot closer...but looking at residual errors one finds that the power1 model may be best overall...will leave that as "exercise for the student" :)
