In this study, a novel algorithm for jointly estimating the linear and quadratic gradients of physical quantities with multiple spacecraft observations based on the least square method has been put forward for the first time. With 10 or more spacecraft constellation measurements as the input, this new algorithm can yield both the linear and quadratic gradients at the barycenter of the constellation. Iterations have been used in the algorithm. The tests on cylindrical flux ropes, dipole magnetic field and modeled geo-magnetospheric field have been carried out. The tests indicate that the linear gradient gained has the second order accuracy, while the quadratic gradient is of the first order accuracy. The test on the modeled geo-magnetospheric field shows that, the more the number of the spacecraft in the constellation, the high the accuracy of the quadratic gradient calculated. However, the accuracy of the linear gradient yielded is independent of the number of the spacecraft. The feasibility, reliability and accuracy of this algorithm have been verified successfully. This algorithm can find wide applications in the design of the future multiple S/C missions as well as in the analysis of multiple point measurement data.