We propose several algorithms for positive polynomial approximation. The main tool is a novel iterative method to compute non negative interpolation polynomials at any order, which is shown to converge under conditions that make it suitable for the numerical approximation of positive functions. Our method is based on the special representations of non negative polynomials provided by the Lukács Theorem, and a key point is the use of Chebyshev polynomials for the initial step of the iterations. Numerical results illustrate the convergence properties of the proposed algorithms, and they are completed with a first application of this technique to the positive discretization of the advection equation.