The approximation of experimental data can be envisaged in the light of fractal interpolation functions defined by iterated function systems. In the particular case of polynomial fractal interpolation functions, the method can be considered as a generalization of the splines of the same kind. That extension is verified under preservation of the smoothness of the function. A bound of the interpolation error by odd degree polynomial fractal interpolation functions is obtained here. The upper limit is also given for high-order derivatives, up to the (2m − 2)th derivative if the polynomials are of degree 2m− 1. The result can be considered a fractal version of Schultz’s theorem for odd degree polynomial splines.