Approximate multidimensional Riemann solvers are essential building blocks in designing globally constraint-preserving finite volume time domain and discontinuous Galerkin time domain schemes for computational electrodynamics (CED). In those schemes, we can achieve high-order temporal accuracy with the help of Runge–Kutta or ADER time-stepping. This paper presents the design of a multidimensional approximate generalized Riemann problem (GRP) solver for the first time. The multidimensional Riemann solver accepts as its inputs the four states surrounding an edge on a structured mesh, and its output consists of a resolved state and its associated fluxes. In contrast, the multidimensional GRP solver accepts as its inputs the four states and their gradients in all directions; its output consists of the resolved state and its corresponding fluxes and the gradients of the resolved state. The gradients can then be used to extend the solution in time. As a result, we achieve second-order temporal accuracy in a single step. In this work, the formulation is optimized for linear hyperbolic systems with stiff, linear source terms because such a formulation will find maximal use in CED. Our formulation produces an overall constraint-preserving time-stepping strategy based on the GRP that is provably L-stable in the presence of stiff source terms. We present several stringent test problems, showing that the multidimensional GRP solver for CED meets its design accuracy and performs stably with optimal time steps. The test problems include cases with high conductivity, showing that the beneficial L-stability is indeed realized in practical applications. © 2023, The Author(s), under exclusive licence to Springer Science+Business Media, LLC, part of Springer Nature.