What is wrong in my Taylor series approximation?

First of all, I apologize for trying to answer this using Python, since I don't have Matlab (even though I know it's a powerful tool).

Your code looks correct to me. The sign in f_curr - b*f_curr_der looks right for approximating $f(x-\delta)$ from values and derivatives at $x$. Your analytical derivative also matches the gradient well. The actual issue is not the sign. It is that the higher-order Taylor terms are not negligible for the high-frequency part of your function. The Taylor expansion

$$f(x-\delta) = \sum_{n=0}^{\infty}\frac{(-\delta)^n}{n!}\,f^{(n)}(x)$$

converges, but the smallness parameter for a Fourier component $\sin(kx)$ is not $\delta$. It is $\delta k$. For your high-frequency term $0.05\sin(60\pi x)$ we have $k=60\pi$ and therefore $\delta\,k = 0.05\cdot 60\pi = 3\pi \approx 9.42$.

You are shifting by roughly $1.5$ wavelengths of the ripple, so no low-order truncation can possibly capture it. The $n$-th term scales like $(\delta k)^n/n!$, which grows until $n\approx\delta k$ and only then starts decaying. In practice we could try $\gtrsim 25$ terms before the truncation error drops below the signal, which is numerically impractical, let's say "a kind of dirty" (high-order derivatives amplify noise).

There are two possible ways to handle this. I use the second one in the code below:

1. Sub-stepping: Apply $M$ smaller shifts $\delta/M$ instead of one large shift $\delta$, so that the relevant smallness parameter becomes $(\delta/M)k$ instead of $\delta k$. Then a low-order Taylor approximation can be applied more safely at each small step.

2. Apply the shift in Fourier space: The shift operator can be written as $e^{-\delta\partial_x}$ and in Fourier space, differentiation becomes multiplication by $ik$, so the shift becomes multiplication by $e^{-ik\delta}$ resulting in $f(x-\delta) = \mathcal{F}^{-1}\!\bigl[e^{-ik\delta}\,\hat f(k)\bigr]$.

Following the latter approach with the following two Jupyter cells (Python), we obtain a good (visual evident) approximation:

import sympy as sp, numpy as np, matplotlib.pyplot as plt
from scipy.fft import fft, ifft, fftfreq

# Define the function symbolically and convert it to a lambda expression
y = sp.symbols('y')
f_sym = sp.cos(2*sp.pi*y) + sp.Rational(1, 20)*sp.sin(60*sp.pi*y)
f = sp.lambdify(y, f_sym, 'numpy')

x = np.linspace(0, 1, 4000, endpoint=False)
b, T = 0.05, 1.0
f_sym

yielding $f\_sym=\sin(60\pi y)/20+\cos(2\pi y)$ as symbolic form of the initial signal ($0.05$ written as $1/20$). And the second Jupyter cell is

k = 2*np.pi * fftfreq(x.size, d=x[1]-x[0])
f_approx = np.real(ifft(np.exp(-1j*k*b*T) * fft(f(x))))
f_exact  = f(x - b*T)

err = np.max(np.abs(f_approx - f_exact))
fig, ax = plt.subplots(figsize=(11, 4))
ax.plot(x, f_exact,  'k',  lw=2.2, label='exact $f(x-bT)$', alpha=0.7)
ax.plot(x, f_approx, 'r--', lw=1.0, label=f'spectral, max|Error|={err:.1e}')
ax.set(xlim=(0.2, 0.7), xlabel='x', ylabel='f(x, t=1)')
ax.legend(); ax.grid(alpha=0.3); plt.tight_layout(); plt.show()

The resulting plot is:

Note that in my Python code fft and ifft denote the discrete Fast Fourier Transform and its inverse (in MATLAB the corresponding functions have the same names fft, ifft).