The research presented in this dissertation is aimed at the development of family of numerical solutions of Two Point Boundary-Value Problems (BVP) in all the branches of engineering sciences with higher precision. In such engineering related boundary-value problems the boundary conditions are specified at two points. The core of this dissertation is focused on the development of a new technique based on quartic non-polynomial spline functions, connecting spline functions values at "mid knots" for the numerical approximations to engineering BVPs and their corresponding values of the fourth-order derivatives. This new approach, which is recently being cited, is so developed that it leads to a family of numerical methods that may be used for computing approximations to the solution of a system of boundary-value problems of third and fourth-order associated with contact and sandwich beam problems. This research work is furthered on the development of a family of higher order numerical solutions to special non-linear third-order boundary-value problems. It is shown that the developed family of higher order methods gives better approximations in contrast to the existing numerical methods. It has further been shown that the existing second, fourth- order and higher order finite-difference and spline functions based methods become special cases of our new technique developed at mid knots presented in this dissertation. Results from the numerical experimentation of renowned problems are also given to illustrate applicability and efficiency of the newly developed family of numerical methods.