We present tests carried out on optical and infrared stellar spectra to evaluate the accuracy of different types of interpolation. Both model atmospheres and continuum normalized fluxes were interpolated. In the first case we used linear interpolation, and in the second linear, cubic spline, cubic-Bezier and quadratic-Bezier methods. We generated 400 ATLAS9 model atmospheres with random values of the atmospheric parameters for these tests, spanning between -2.5 and +0.5 in [Fe/H], from 4500 to 6250 K in effective temperature, and 1.5 to 4.5 dex in surface gravity. Synthesized spectra were created from these model atmospheres, and compared with spectra derived by interpolation. We found that the most accurate interpolation algorithm among those considered in flux space is cubic-Bezier, closely followed by quadratic-Bezier and cubic splines. Linear interpolation of model atmospheres results in errors about a factor of two larger than linear interpolation of fluxes, and about a factor of four larger than high order flux interpolations.