An interpolating spline which interpolates positive function values is not necessarily positive for all arguments; a typical example of this occurs in the spline interpolation of the normal density function of probability theory. We describe an algorithm which produces non-negative interpolating splines provided the function to be interpolated is positive. We start with the interpolating natural spline and redefine it in those intervals for which the spline has negative values; we do this by adding few extra knots to those intervals. First we show that for any choice of extra knots it is possible to construct a spline which is non-negative. We call such a spline feasible. Second we show that within the set of feasible splines there is a best spline, in the sense that the functional which reflects the strain energy is minimized. Finally, by using a necessary optimality criterion we obtain an optimal spline which minimizes the strain energy also with respect to the free knots. The resulting formulae give explicit results and are very simple. Our numerical results show that the final step of selecting the interior knots appropriately gives very satisfactory results in a very efficient manner. The algorithm can be applied also to produce splines which stay below upper or stay above lower constant bounds. It is also very useful as a starting procedure for finding the globally best spline.