Abstract
We use quadratic spline functions to develop a numerical method for computing approximate solution of a system of second order boundary value problems associated with obstacle, unilateral and contact problems. The present method outperforms other available techniques and it has extra added advantage of approximating the first derivative of the solution. Numerical example is presented to illustrate the applicability of the new method.