A computational method for the solution of differential equations is proposed. With this method an accurate approximation is built by incremental additions of optimal local basis functions. The parallel direct search software package (PDS), that supports parallel objective function evaluations, is used to solve the associated optimization problem efficiently. The advantage of the method is that, although it resembles adaptive methods in computational mechanics, an a priori grid is not necessary. Moreover, the traditional matrix construction and evaluations are avoided. Computational cost is reduced while efficiency is enhanced by the low-dimensional parallel-executed optimization and parallel function evaluations. In addition, the method should be applicable to a broad class of interpolation functions. Results and global convergence rates obtained for one- and two-dimensional boundary value problems are satisfactorily compared to those obtained by the conventional Galerkin finite element method