This paper introduces methods and a novel toolbox that efficiently integrates any high-dimensional Neural Mass Models (NMMs) specified by two essential components. The first is the set of nonlinear Random Differential Equations of the dynamics of each neural mass. The second is the highly sparse three-dimensional Connectome Tensor (CT) that encodes the strength of the connections and the delays of information transfer along the axons of each connection. Semi-analytical integration of the RDE is done with the Local Linearization scheme for each neural mass model, which is the only scheme guaranteeing dynamical fidelity to the original continuous-time nonlinear dynamic. It also seamlessly allows modeling distributed delays CT with any level of complexity or realism, as shown by the Moore-Penrose diagram of the algorithm. This ease of implementation includes models with distributed-delay CTs. We achieve high computational efficiency by using a tensor representation of the model that leverages semi-analytic expressions to integrate the Random Differential Equations (RDEs) underlying the NMM. We discretized the state equation with Local Linearization via an algebraic formulation. This approach increases numerical integration speed and efficiency, a crucial aspect of large-scale NMM simulations. To illustrate the usefulness of the toolbox, we simulate both a single Zetterberg-Jansen-Rit (ZJR) cortical column and an interconnected population of such columns. These examples illustrate the consequence of modifying the CT in these models, especially by introducing distributed delays. We provide an open-source Matlab live script for the toolbox.