How to find ePSF for multiple dithered undersample images?

[Eririf]

I want to use ePSF to refine star center positions, then feed to Scamp for better wcs to do SWarp. But im not clear how EPSFBuilder work on dither images, LinkedEPSFStar seems work but details not provided in tutorial.

fitsfiles=glob.glob(workload+'*.fits')
from astropy.visualization import simple_norm
i=0
fitsfile=fitsfiles[i]
hdu=fits.open(fitsfile)
data=hdu[0].data
norm = simple_norm(data, 'sqrt', percent=99.)
plt.imshow(data, norm=norm, origin='lower', cmap='viridis')
hdu.close()
from photutils.detection import find_peaks
peaks_tbl = find_peaks(data, threshold=7000.)  
peaks_tbl['peak_value'].info.format = '%.8g'  # for consistent table output  
peaks_tbl=peaks_tbl[(peaks_tbl['peak_value']<50000)]

size = 5
hsize = (size - 1) / 2
x = peaks_tbl['x_peak']  
y = peaks_tbl['y_peak']  
mask = ((x > hsize) & (x < (data.shape[1] -1 - hsize)) &
        (y > hsize) & (y < (data.shape[0] -1 - hsize)))  

from astropy.table import Table
stars_tbl = Table()
stars_tbl['x'] = x[mask]  
stars_tbl['y'] = y[mask]  

from astropy.stats import sigma_clipped_stats
mean_val, median_val, std_val = sigma_clipped_stats(data, sigma=2.)  
data -= median_val  

from astropy.nddata import NDData
nddata = NDData(data=data)  

from photutils.psf import extract_stars
stars = extract_stars(nddata, stars_tbl, size=(11,11))  

from photutils.psf import EPSFBuilder
epsf_builder = EPSFBuilder(oversampling=2, maxiters=5,smoothing_kernel='quartic',
                           progress_bar=False)  
epsf, fitted_stars = epsf_builder(stars)  
import matplotlib.pyplot as plt
from astropy.visualization import simple_norm
norm = simple_norm(epsf.data, 'log', percent=99.)
plt.imshow(epsf.data, norm=norm, origin='lower', cmap='viridis')
plt.colorbar()