We address in this paper a new approach for fitting spatiotemporal models with application in disease mapping using the interaction types I, II, III, and IV proposed by [1]. When we account for the spatiotemporal interactions in disease-mapping models, inference becomes more useful in revealing unknown patterns in the data. However, when the number of locations and/or the number of time points is large, the inference gets computationally challenging due to the high number of required constraints necessary for inference, and this holds for various inference architectures including Markov chain Monte Carlo (MCMC) [2] and Integrated Nested Laplace Approximations (INLA) [3]. We re-formulate INLA approach based on dense matrices to fit the intrinsic spatiotemporal models with the four interaction types and account for the sum-to-zero constraints, and discuss how the new approach can be implemented in a high-performance computing framework. The computing time using the new approach does not depend on the number of constraints and can reach a 40-fold faster speed compared to INLA in realistic scenarios. This approach is verified by a simulation study and a real data application, and it is implemented in the R package INLAPLUS and the Python header function: inla1234().