In this thesis, I discuss the work I have contributed towards the development andapplication of CSbosonsCpp, a software package written in C++11 and CUDA, designed to simulate systems of many interacting bosons at finite temperature. The simulation method uses Feynman’s path integral in the basis of coherent states fields as the theoretical underpinning and complex Langevin sampling to calculate relevant thermodynamic quantities, such as number of particles, chemical potential, and free energy, of systems of many interacting bosons. First, I review the theory and numerical methods upon which this software is built. I briefly compare the performance of this software as applied to a system of many interacting bosons with a GPU-accelerated Path Integral Monte Carlo (PIMC) software package called cudapimc maintained by Prof. Adrian Del Maestro at the University of Tennessee, Knoxville and Benjamin Himberg at the University of Vermont. PIMC is the current premier method for simulating many-body quantum systems at finite temperature. I finally detail two studies in which we apply CSbosonsCpp to systems for which which PIMC is not a practical choice. The first system is cloud of interacting bosons mechanically rotated in the laboratory frame. PIMC is a non-ideal choice for simulating this system due to a sampling issue known as the sign problem. The second is a thermodynamic engine that uses a Bose-Einstein condensate of 7Li as the working fluid. The system size of half a million particles is untenable for current PIMC algorithms.