To have all the samples in the interpolated waveform land on the original samples implies using an integer interpolation rate (otherwise time domain distortion would need to be introduced, and I can't think of an application where that would make sense to do).
For integer interpolation, a very simple approach is to use the DFT with zero padding since the samples chosen will be exact, but this will have more deviation in between then other optimized interpolation methods further detailed below. This approach is very simple, so if the deviation between samples is of no concern then this approach would be favored and would be as follows:
Compute the FFT of the samples without zero padding.
Scale the result by $N$, where $N$ is the total number of original samples, zero pad the result and then compute the inverse FFT. (Due to reciprocity the fft and ifft operations can be swapped).
In MATLAB or Octave this is:
N = length(test);
interp = 5;
result = ifft(N*fft(sequence), interp*N);
Other optimized interpolation methods using zero-insert and subsequent image filtering would also preserve the values, and with post processing the image filtering can be done with a zero-phase filter using
filtfilt() and for an optimized solution in the least squares sense, together with multiband filter coefficients using
firls(). The zero-insert and filtering approach is further detailed here: What is the impulse response used in an interpolation filter when upsampling?
If your samples are non-uniform to begin with, first resample to a uniform rate, where this post already addresses that approach: What is an algorithm to re-sample from a variable rate to a fixed rate?
But if the new samples do not fall on the same time locations (which will occur if there isn't an integer relationship in the duration between the original samples) then the interpolated samples should not be expected to match the original samples when they occur at new time samples that were not in the original sequence.