Simulating turbulence in the tokamak boundary, which is crucial for understanding and optimizing future fusion power plants performance, is computationally challenging due to the complex, multi-physics, and multi-scale nature of plasma dynamics. In the boundary region, plasma interacts with neutral particles recycled from the walls, leading to nonlinear phenomena that control confinement and heat exhaustion. Accurately capturing these coupled processes requires solving large-scale, stiff, and non-local problems. This thesis contributes to the development of efficient numerical algorithms based on low-rank linear algebra to accelerate plasmaâ neutral simulations. The key idea is to exploit the low-dimensional structure underlying the data and operators, which can often be approximated in reduced subspaces, thereby lowering memory and computational cost. The proposed methods are implemented in GBS, a three-dimensional, flux-driven turbulence code for self-consistent plasma and kinetic neutral modeling in the tokamak boundary. GBS solves the drift-reduced Braginskii equations coupled with a kinetic neutral model that is solved deterministically through integration along characteristics. The first contribution of this thesis introduces a subspace acceleration method for expediting the solution of sequences of large-scale linear systems, such as the ones arising from the numerical discretization of time-dependent partial differential equations coupled with algebraic constraints. By leveraging the history of previous solutions, accurate initial guesses for an iterative solver are constructed through reduced-order projection methods combined with randomized linear algebra techniques, drastically reducing the number of iterations needed for convergence. Theoretical analysis highlights that smooth temporal evolution ensures high approximation quality, while large-scale simulations of the plasma boundary demonstrate significant speed-up, through faster solution of the resultant linear systems. The second contribution concerns the acceleration of the solver used to simulate the kinetic neutral dynamics. When the neutral dynamics are described deterministically via the method of characteristics, the resulting integral formulation leads to dense linear systems which are too costly to assemble and solve in practice. This thesis employs hierarchical matrix approximations to represent the corresponding operators in a data-sparse format, leading to a computational complexity nearly linear in the number of unknowns, in contrast with a quadratic scaling of the dense approach. The hierarchical matrix solver achieves over 90% reduction in computation time and memory load, enabling high-resolution simulations of neutral dynamics consistent with the plasma grid. These developments advance the efficiency and accuracy of tokamak boundary simulations, enabling detailed and computationally tractable modeling of plasmaâ neutral interactions, which are essential for the design and operation of future fusion devices.