An inverse method is introduced to construct benchmark problems for the numerical solution of initial value problems. Benchmark problems constructed through this method have a known exact solution, even though analytical solutions are generally not obtainable. The solution is constructed such that it lies near a given approximate numerical solution, and therefore the special case solution can be generated in a versatile and physically meaningful fashion and can serve as a benchmark problem to validate approximate solution methods. A smooth interpolation of the approximate solution is forced to exactly satisfy the differential equation by analytically deriving a small forcing function to absorb all of the errors in the interpolated approximate solution. A multi-variable orthogonal function expansion method and computer symbol manipulation are successfully used for this process. Using this special case exact solution, it is possible to directly investigate the relationship between global errors of a candidate numerical solution process and the associated tuning parameters for a given code and a given problem. Under the assumption that the original differential equation is well-posed with respect to the small perturbations, we thereby obtain valuable information about the optimal choice of the tuning parameters and the achievable accuracy of the numerical solution. Illustrative examples show the utility of this method not only for the ordinary differential equations (ODEs) but for the partial differential equations (PDEs).