Quartic non-polynomial spline functions are used to develop numerical methods for computing approximations to the solution of systems of second-order boundary-value problems associated with obstacle, unilateral, and contact problems. It is shown that the present methods give approximations, which are better than those produced by other collocation, finite-difference and spline methods. Convergence analysis of the method for obstacle problems is discussed. Numeric example's given to illustrate practical usefulness of the new approach.