Your problem is, you are trying to use LEAST squares. But look carefully at your data.
xdata = [0.02 0.05 0.1 0.2 0.3 0.4 0.5];
ydata = [1.7935e-006 5.6359e-006 10.8276e-006 38.1193e-006 60.4152e-006 65.8997e-006 62.0615e-006];
Do you see that some of your data is nearly 2 powers of 10 larger then the rest? You plot it on a log scale, so you do't really think of it. It looks like just a straight line on those axes.
So what happens is error in some of your data points are WAY more important than are errors in the others. This is a case where you do NOT want to use a simple least squares to estimate those coefficients. You might decide to use a weighted least squares. A simple alternative is to use a straight line fit, to the log of your data.
Why is that a good idea? Because in the case you have, you really do not want to use an error structure that is additive, with a uniform variance on all data points. Instead, you want to use a proportional error structure, s multiplicative errors. What happens when you take the log? It turns what you have, thus a lognormally distributed error, into normally distributed, additive error.
xlog = log(xdata);
ylog = log(ydata);
P1 = fit(xlog',ylog','poly1')
plot(P1,xlog,ylog,'o')
I used fit there, because it makes the plotting really easy.
But you should see how well it worked here, giving you the model you were hoping to find.
Now we need to transform the model back into a and b.
b = P1.p1
a = exp(P1.p2)
The distinction I pointed out here is an important one.
