def METHOD_NAME(im, psfmodel, cyx, scaling='peak'):
    """subtract a PSF model from the image
    Parameters
    ----------
    im : array, mandatory
        image
    psfmodel : array, mandatory
        PSF model image
    cyx : tuple, mandatory
        target centroid coordinates `(y, x)`
    scaling : str, optional, default: `peak`
        select scaling method (`peak` fits the height of PSF model to
        the target's image)
    Examples
    --------
    >>> from astropy.io import fits
    >>> from sbpy.imageanalysis import centroid, PSFSubtraction # doctest: +SKIP
    >>> import matplotlib.pyplot as plt # doctest: +SKIP
    >>> hdu = fits.open('test.fits') # doctest: +SKIP
    >>> cyx = centroid(hdu[0].data, (21, 45)) # doctest: +SKIP
    >>> psfmodel = PSFSubtraction.create_psfmodel(hdu[0].data, cyx) # doctest: +SKIP
    >>> diff = PSFSubtraction.subtract_psfmodel(hdu[0].data, psfmodel, cyx) # doctest: +SKIP
    >>> plt.imshow(diff) # doctest: +SKIP
    not yet implemented
    """