The computation of the reachable set of states of a given dynamic system is an important step to verify its safety during operation. There are different methods of computing reachable sets, namely interval integration, capture basin, methods involving the minimum time to reach function, and level set methods. This work deals with interval integration to compute subpavings to over or under approximate reachable sets of low dimensional systems. The main advantage of this method is that, compared to guaranteed integration, it allows to control the amount of over-estimation at the cost of increased computational effort. An algorithm to over and under estimate sets through subpavings, which potentially reduces the computational load when the test function or the contractor is computationally heavy, is implemented and tested. This algorithm is used to compute inner and outer approximations of reachable sets. The test function and the contractors used in this work to obtain the subpavings involve guaranteed integration, provided either by the Euler method or by another guaranteed integration method. The methods developed were applied to compute inner and outer approximations of reachable sets for the double integrator example. From the results it was observed that using contractors instead of test functions yields much tighter results. It was also confirmed that for a given minimum box size there is an optimum time step such that with a greater or smaller time step worse results are obtained.