If you're working with hg38, you could do something like the following to get counts of reads from a BAM file in 50kb bins (requires UCSC Kent tools and BEDOPS):

$ fetchChromSizes hg38 | awk -v OFS="\t" '($0!~/.*_.*/){print $1, "0", $2; }' | sort-bed - > hg38.bed
$ bedmap --echo --count --delim '\t' <(bedops --chop 50000 hg38.bed) <(bam2bed < reads.bam) > hg38.50kbReadDensity.bed

You could adjust the width of this window by changing 50000 to some other value.

If you have two BED files, one of which that contains signal in the score column, you could instead do this:

$ bedmap --echo --count --delim '\t' <(bedops --chop 50000 hg38.bed) reads.bed > hg38.50kbReadDensity.bed

Either way, you could then split this file by chromosome:

$ for chr in `bedextract --list-chr hg38.50kbReadDensity.bed`; do echo ${chr}; bedextract ${chr} hg38.50kbReadDensity.bed > hg38.50kbReadDensity.${chr}.bed; done

Then plot the signal in each per-chromosome file to a per-chromosome PDF:

$ for chr in `bedextract --list-chr hg38.50kbReadDensity.bed`; do echo ${chr}; ./plotDensity.Rscript hg38.50kbReadDensity.${chr}.bed hg38.50kbReadDensity.${chr}.pdf "50kb Read Density - ${chr}"; done

The plotDensity.Rscript script could look like this:

#!/usr/bin/env Rscript
args = commandArgs(trailingOnly=TRUE)
inFn = args[1]
outFn = args[2]
title = args[3]
d <- read.table(inFn, col.names=c("chr", "start", "stop", "density"), sep="\t", colClasses=c("character", "numeric", "numeric", "numeric"))
xRange <- range(0, d$stop) 
yRange <- range(0, d$density)
pdf(outFn)
plot(xRange, yRange, type="n", xlab="Position", ylab="Density", main=title)
lines(d$start, d$density)
dev.off()