function a = mysvdfit(x, y)
%MYSVDFIT SVD least squares fit of (x,y) data.
%          Return the coefficients of the polynomial
%          p(x) = a(1) + a(2)*x + a(3)*x^2 + ... + a(n+1)*x^n
    np = length(x);
    A = [ones(1, np)' x' x'.^2];
    [U, S, V] = mysvd(A);
    r = nnz(diag(S));
    Uhat = U(:,1:r);
    Shat = S(1:r,1:r);
    z = Shat\Uhat'*y';
    a = V*z;
end

