This paper presents a new method for polynomial approximation using the fusion of value and derivative information emanating from different sources, i.e., sensors. Therefore, the least-squares error in both domains is simultaneously minimized. A covariance weighting is used to introduce a metric between the value and derivative domain, to handle different noise behaviour. Based on a recurrence relation with full re-orthogonalization, a weighted polynomial basis function set is generated. This basis is numerically more stable compared to other algorithms, making it suitable for the approximation of data with high degree polynomials. With the new method, the fitting problem can be solved using inner products instead of matrix-inverses, yielding a computational more efficient method, e.g., for realtime approximation.A Monte Carlo simulation is performed on synthetic data, demonstrating the validity of the method. Additionally, various tests on the basis function set are presented, showing the improvement on the numerical stability.