The HBC Training link gives most details you'll need. The basic steps are:
- obtain normalized coverage files of your CUT&RUN samples (bigWig files)
- use these bigWig files together with the MACS BED files to generate the matrices needed for heatmaps and coverage plots.
If you want to force a specific order of the regions in the BED file (e.g. from "highest" to "lowest"), you'll have to turn off the clustering for the heatmap. Or you can sort the regions by the mean expression values for either sample (as computed by deepTools rather than MACS2); check out the optional arguments for plotHeatmap: https://deeptools.readthedocs.io/en/develop/content/tools/plotHeatmap.html#Optional%20arguments