We present a computational strategy for the evaluation of multidimensional integrals on hyper-rectangles based on Markovian stochastic exploration of the integration domain while the integrand is being morphed by starting from an initial appropriate profile. Thanks to an abstract reformulation of Jarzynski’s equality applied in stochastic thermodynamics to evaluate the free-energy profiles along selected reaction coordinates via non-equilibrium transformations, it is possible to cast the original integral into the exponential average of the distribution of the pseudo-work (that we may term “computational work”) involved in doing the function morphing, which is straightforwardly solved. Several tests illustrate the basic implementation of the idea, and show its performance in terms of computational time, accuracy and precision. The formulation for integrand functions with zeros and possible sign changes is also presented. Program summaryProgram title: JEMDICatalogue identifier: AEWO_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AEWO_v1_0.htmlProgram obtainable from: CPC Program Library, Queen’s University, Belfast, N. IrelandLicensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.htmlNo. of lines in distributed program, including test data, etc.: 81666No. of bytes in distributed program, including test data, etc.: 1767836Distribution format: tar.gzProgramming language: C++.Computer: PC and Macintosh.Operating system: Unix.RAM: 1 024 000 000 bytesClassification: 4.11, 4.13.External routines: SPRNG included in the distribution file.Nature of problem:Integration of general functions of many variablesSolution method:JEMDI is a library that provides methods for the evaluation of general multidimensional integrals on hyper-rectangles based on Markovian stochastic exploration of the integration domain while the integrand is being morphed by starting from an initial flat profile. Employing a revised abstract reformulation of Jarzynski’s equality, the original integral is cast into the exponential average of the distribution of the pseudo-work involved in doing the morphing of the integrand function.Restrictions:The method is efficiently applicable to the integration of functions of a few hundreds of variables.Running time:Extremely variable, depending on the complexity of the function to integrate and the number of variables. The examples included in the distribution take about 3 min to run.