Support PSFEx variable psf

Hi,

It would be nice if support for PSFEx psf file can be added in photutils (at least in limited capacity).

The structure of the variable psf file is a fits files containing two extensions and the structure is as follows -
The 2nd extension (hdul[1]) contains the relevant header and PSF data.

The header contains the following import parameters -
1. PSFvariation_dependent variable (say x, and y)
2. variation degree
3. Oversampling factor (float)

The data is an array shaped (1, 1, n, psf_x, psf_y) where n is the coefficient for the polynomials.
e.g. for a 3rd order variation, n=10.
The coefficients, according to the variation degree are in the following order -
vardeg=0 -> [1]
vardeg=1 -> [1, x, y]
vardeg=2 -> [1, x, x**2, y, y*x, y**2]
vardeg=3 -> [1, x, x**2, x**3, y, y*x, y*x**2, y**2, y*x**2, y**3]
and so on.

The variation can depend on any image parameter (e.g. flux) or any other fits parameter (e.g. temperature, focus) but that might be difficult to implement. Therefore, I am suggesting a limited support that just covers x and y coordinates of the image.

I looked around the code, and It could just be a subclass of the FitttableImageModel. We would need to implement a new `evaluate' method.
For every x_0 and y_0 we can calculate a coefficient array -> multiply it with psf_array -> sum all arrays together to get the final PSF.
However, I am unsure how it will affect the other internal function of the class.

Start by reading the FitttableImageModel implementation and its evaluate method, then inspect how PSFEx FITS extensions, headers, and data arrays are represented. Define the x/y-only variable-PSF behavior from the stated coefficient ordering and shape; done means the extension data and oversampling information produce the expected PSF across image coordinates.

Written by the indexing model from the issue text.