Simulating a multipath channel
The upsampling solution adds additional difficulties, as mentioned in some posts here. Please consider this solution..

The Discrete Fourier Transform (DFT) properties might support you. Take a look at the time shifting property. It states that $$\mathcal{F}\{x[k-\tau]\} = e^{-j\omega\tau}\mathcal{F}\{x\}$$ I implemented a possible use of the property in your problem in the following octave script. A picture with the result is also attached. A cos signal with 3[Hz] is sampled at 10[Hz], and shifted 33[ms] with the property. I hope it helps.

% sampled signal:
fs          = 10;
Ts          = 1/fs;
T           = 0:Ts:1;
nT          = length(T);
sig_sampled = 0.2 * cos(2*pi*3*T);

% just to plot, signal with high sample rate:
Tb          = 0:(1/1000):1;
sig_better  = 0.2 * cos(2*pi*3*Tb);

% delayed signal where tau is not a multiple of Ts:
% (using Fourier properties):
tau              = 0.033;
w                = 2*pi*(((-nT/2):((nT/2)-1))/nT)*fs;
fft_sig          = fft(sig_sampled);
fft_shift        = exp(-sqrt(-1)*tau*ifftshift(w)) .* fft_sig;
sig_delayed      = ifft(fft_shift);
sig_better_shift = 0.2 * cos(2*pi*3*(Tb-tau));

figure(1); clf;
plot(
    Tb, sig_better, 'k-',
    T,sig_sampled, 'k-x',
    Tb, sig_better_shift, 'r-',
    T, sig_delayed, 'r-x'
);
set(gca, 'ylim', 0.3*[-1,1])
grid on
legend('original signal', 'sampled', 'shifted signal', 'shifted sampled signal')
print('-djpg', '-r 200', 'fig.jpg')