I would use either the envelope function (also ‘abs(hilbert(signal))’) or the Savitzky-Golay filter function sgolayfilt.
To use either of those, you will probably have to do it twice, first on the signal as it currently exists, then on the negative of the signal. Use the ‘negative of the negative’ result to plot the lower boundary. You can add or subtract small offsets to get the result you want.
Everything in signal processing requires experimentation, so you will most likely have to tweak the function parameters a few times.