Of the methods you've already tried, the "looking for peaks in the spectrum" to me sounds like the most promising, but:
My guess is that your "spectrum" is actually a PSD estimate done with the FFT. Which is a fine method of getting an understanding of the spectrum – but as you probably noticed, the "peaks" will not be very clear.
Think about it: The whole principle of FM is to shift a carrier's frequency around. So, looking for energy in a single FFT bin won't work overly well, because these bins will be too narrow to "contain" all the energy of a single station (unless you choose very low FFT length -> wide bins, but then you can't tell stations apart).
A different way of looking at the FFT (which is just an implementation of the DFT) than being a (base) transform between time and frequency domain is understanding it as filter bank of sinc-shaped filters. So, if your radio stations' spectra happened to have sinc shape, and be $\frac{f_\text{sample}}{N_\text{FFT}}$ spaced, this would work indeed very well.
Sadly, the description doesn't fit FM radio stations. But: with a less specific filter bank, we might be able to find channels!
First observation is that if it is somewhat wide in spectrum, and in the UHF band between 87 and 108 MHz, it's an FM station, or illegal. So, let's not waste too much thought on the specific spectral shape of an FM station¹; instead, let's just
- use a filter that lets through about 160 kHz of bandwidth
- Calculate the power of all that signal that passes the filter (simply by magnitude-squaring the samples and averaging the output)
- Take that filter+detector and apply it to all potential channels
Filtering Now, the design of that filter (1.) might be rather simplistic: I'll take an arbitrary real-valued low-pass filter with a cutoff frequency of 80 kHz; that will let through the 160 kHz between -80 kHz and +80 kHz around its center frequency.
Magnitude squaring Remember that digital samples are proportional to the voltage at the ADC, but will be digital here. |·|² is simply a way of calculating a number that is proportional to the power of the signal – and since we're in no way interested in absolute powers, but only in what is "stronger" than some arbitrary threshold, this will do.
Applying the filter all over the spectrum The old-school spectrum analyzer method would be to take a mixer, tune it sequentially to different frequencies, apply the filter to the downmixed signal, and calculate the power. In fact, we can do something like that! But then, we only get info about a single station at a time. That might or might not be good enough for us, and it's certainly as good as the average car radio does (which also scans!).
Now, we're proud applicants of DSP, so we might as well go back to the filter bank idea: Let's just take that filter, replicate it, shifted by 50 kHZ (which I presume is the FM station frequency raster), and apply those 420 filters simultaneously!
Well, that would require that you "see" the full 21 MHz of FM-allocated spectrum at once, and your RTL-dongle won't do that. Instead, let's settle with the 2 MHz that your dongle are able to sample at once. That's 40 filters. I'll come back later and have a signal processing flow graph and code for you.
¹ Bessel functions, by the way. You don't want to analytically deal with those.
² $\frac{(108 - 87)\,\text{MHz}}{50\,\text{kHz}}=\frac{21\cdot10^3}{5\cdot10^1}=420$
