Context. Electromagnetic waves arise in many areas of physics. Solutions are difficult to find in the general case. Aims. We numerically integrate Maxwell equations in a 3D spherical polar coordinate system. Methods. Straightforward finite difference methods would lead to a coordinate singularity along the polar axis. Spectral methods are better suited for such artificial singularities that are related to the choice of a coordinate system. When the radiating object rotates like a star, for example, special classes of solutions to Maxwell equations are worthwhile to study, such as quasi-stationary regimes. Moreover, in high-energy astrophysics, strong gravitational and magnetic fields are present especially around rotating neutron stars. Results. To study such systems, we designed an algorithm to solve the time-dependent Maxwell equations in spherical polar coordinates including general relativity and quantum electrodynamical corrections to leading order. As a diagnostic, we computed the spin-down luminosity expected for these stars and compared it to the classical or non-relativistic and non-quantum mechanical results. Conclusions. Quantum electrodynamics leads to an irrelevant change in the spin-down luminosity even for a magnetic field of about the critical value of 4.4 × 10 9 T. Therefore the braking index remains close to its value for a point dipole in vacuum, namely n = 3. The same conclusion holds for a general-relativistic quantum electrodynamically corrected force-free magnetosphere.