Estimate phase using STFT in Matlab?
Jumping between bins of STFTs will lead to the discontinuities seen in your phase estimates.  Interpolating both frequency and phase from a sequence of STFTs is possible, but requires non-standard techniques, such as fftshifts and 2D interpolators.

If you know the frequency range of your desired slow oscillation, then a bandpass filter in conjunction with a PLL follower might provide a better looking phase result using more standard DSP techniques.