How to remove spikes in solution and produce smooth interpolation with scipy?
I don't work with scipy, but from what I've gathered, some things stood out to me, and could possibly be what's causing problems.

Your call to plt.show() which displays the data happens before you filter out the outliers with medfilt(), so the corrected data might not appear in your plot
The median filter from what I gather doesn't remove outliers from your data, instead it resets each data point with the median value amongst its k-neighbors.
With this in mind, I have two suggestions, (1) your median filter window might be too small, and that is causing the outliers to not be removed. Try setting it yourself using medfilt(self.n_, k_size=5), it defaults to 3 so try odd numbers larger than 3. (2) Given that you're not losing data points from using the medfilter, you might not need the lines that follow it which try to interpolate data that was presumably removed.
I'm writing a numerical solution to a partial integro-differential equation, and I need it to run quickly so I found that scipy.integrate.simps is best, but it's not always 100% accurate and produces the spikes in [1]. My solution was to remove them with scipy.signal.medfilt and then interpolate over the gaps with an interpolator (I've tried CubicSpline, PChipInterpolator, scipy.interp1d, akima,...) but all of them produce little "hiccups" in the solution that can be seen at y=0.1, (produced with 3rd order butterworth filter) and these errors grow as the solution is evolved. How do I remove the spikes and get a simple, smooth interpolation over the gaps? Thanks!