Conventional spline procedures have proven to be effective and useful for estimating smooth functions. However, these procedures find piecewise and inhomogeneous smooth functions difficult to handle. Conventional spline procedures often yield an overly-smooth curve in regions where the true function or its derivative are discontinuous. To fully realize the potential of the spline methodology, this article proposes a function estimation procedure that uses adaptive free-knot splines and allows for multiple knots to be replaced at the same location and between design points. The article also proposes an associated algorithm for implementation in non-parametric regression. The proposed knot selection scheme uses a data-adaptive model selection criterion and an evolutionary algorithm that incorporates certain key features of simulated annealing. The evolutionary algorithm accurately locates the optimal knots, while the data-driven penalty guards against selection errors when searching through a large candidate knot space. The algorithm stochastically yields the globally optimal knots. The simulations suggest that the procedure performs competitively well against alternative methods and has a substantial advantage in relation to non-smooth and piecewise smooth functions.