In this paper, a method to exactly sample the trajectories of inverse subordinators (in the sense of the finite-dimensional distributions), jointly with the undershooting or overshooting process, is provided. The method applies to general strictly increasing subordinators. The (random) running times of these algorithms have finite moments and explicit bounds for the expectations are provided. Additionally, the Monte Carlo approximation of functionals of subdiffusive processes (in the form of time-changed Feller processes) is considered where a central limit theorem and the Berry-Esseen bounds are proved. The approximation of time-changed Itô diffusions is also studied. The strong error, as a function of the time step, is explicitly evaluated demonstrating the strong convergence, and the algorithm's complexity is provided. The Monte Carlo approximation of functionals and its properties for the approximate method is studied as well. An application of our algorithms in the context of weak ergodicity breaking of subdiffusion is also discussed.