It looks to me as if you are defining an arc length. However, for an arc length, f0 would have to be the derivative of a function. Which function?
f0 = ((q(1)*a^5)+(q(2)*a^4)+(q(3)*a^3)+(q(4)*a^2)+(q(5)*a)+q(6));
It is not the derivative of the function expressed through the q coefficients, at least not at that stage. Let's look further back
well, that certainly looks like it might be a derivative. But is it?
No! The output of polyfit() is numeric, and diff() applied to a numeric vector is the finite difference operator, not the derivative. You are computing
q = f1(2:end) - f1(1:end-1)
not taking the derivative of a polynomial implied by f1.
I suspect that you should be doing
g0 = sqrt(1 + diff(f0,a)^2);
You are unlikely to get a closed form solution. Closed form solutions might be possible in terms of the Elliptic Integral if diff(f0,a)^2 was degree no more than 4 (which would require that diff(f0,a) be of degree no more than 2, which would be for the case of fitting a cubic)