Nonpolynomial Spline Methods (NPSMs) are presented for solving single and system of boundary value problems (BVPs) of orders two, four, six and an 2mth order BVPs respectively. Also, the numerical solution of singular second order BVPs is considered. Optimal order approximations to the solutions of these boundary value problems are obtained. The matrix properties of these BVPs, arising from the discretization of systems of BVPs by NPSMs, are discussed. Sufficient conditions for the optimal convergence for NPSMs are obtained. Also, NPSMs generalize other recent finite difference methods and polynomial spline methods through arbitrary choices of some defined parameters in these methods. Numerical results demonstrated by the NPSMs are shown to produce smaller approximation errors and the superiority of our methods over other recent methods such as finite difference methods, polynomial spline methods and spline collocation methods. These new proposed methods are very encouraging in dealing with boundary value problems.