The problem of cubic spline interpolation of functions with large gradients in the boundary layer is considered. It is assumed that the decomposition in the form of the sum of the regular and singular components is valid for the interpolated function. This decomposition is valid for the solution of a boundary value problem for a second order ordinary differential equation with a small parameter e at the highest derivative. An overview of approaches to the construction of splines, the error of which is uniform with respect to the small parameter e, is given. The approaches are based on the use of Shishkin and Bakhvalov meshes, which are thickened in the boundary layer region. An alternative approach based on a modification of a cubic spline with fitting to the singular component of the interpolated function is also considered. The comparison of the accuracy of the applied approaches with the performance of computational experiments is carried out.