Context. Photometric surveys require precise point spread function (PSF) characterization, as the PSF varies across filters and plays a crucial role in achieving accurate photometry, particularly in low surface brightness (LSB) studies. Such studies include the analysis of faint and extended structures like galaxy halos, tidal features, diffuse circumgalactic emission, and intracluster light. However, the small PSF models produced by default data-reduction pipelines are optimized mainly for compact or barely resolved sources, which makes it challenging to analyze regions near bright stars, often rendering those areas unusable. Aims. We aim to demonstrate the feasibility of subtracting the extended PSF from each J-PLUS DR3 exposure prior to sky subtraction, an approach that has not yet been explored in wide surveys, and enable a comprehensive analysis of its behavior across detector position and time. Methods. To build an extended non-parametric PSF, we selected three different ranges of stars to create the central, middle, and outer regions in exposures within ~2.5 hours of the target. These components were then combined to generate a final PSF for each exposure and filter, spanning 15 mag arcsec−2 in surface brightness and 4 arcmin in radius in the broad bands. Results. In narrowband filters, the J-PLUS PSF exhibits two rings, whereas in broadband filters, only one ring is observed. Additionally, the position of the ring shifts with the filter wavelength in the following manner: as the filters become redder, the ring radius increases. We find that the precision of sky subtraction can be greatly improved with a PSF-subracted image, and out to a radius of 4 arcmin, there is no significant variation in the extended PSF observed as a function of time or position in the field of view. The radial profile of NGC 4212 (which is close to a star) is also studied before and after PSF subtraction as a demonstration of the effect. We developed a novel method to determine the central coordinates of saturated stars, and we classified stars without using Gaia magnitudes. Additionally, mirror reflections were automatically detected and masked. Furthermore, in combining different stars and various components of the PSF, we avoided the use of a fixed radius by introducing a new method that does not depend on radial measurements. Conclusions. Accurate characterization of the extended PSF and its subtraction improves sky subtraction, increases the effective area of the survey by about 10%, and enables the study of extended large LSB features in wide area surveys such as J-PLUS. Our pipeline is published as free software (GNU GPLv3), and it can be customized to suit other surveys such as J-PAS, where its impact will be even greater due to its depth. This paper is fully reproducible and produced from Commit 21523b8.