Assuming $f(x)=e^{-x^2}$ for $x$ in $[-10,10]$, I have tried the following:


1. Fourier transform $\mathcal{F}$: $\frac{d}{dx}$ can be diagonalized as $\mathcal{F}^{-1} i\omega \mathcal{F}$. Therefore, $\sin(\frac{d}{dx}) f(x) = \sin(\mathcal{F}^{-1} i\omega \mathcal{F})f(x)=\mathcal{F}^{-1}\sin(i\omega)\*\mathcal{F}f(x)=\mathcal{F}^{-1}\sinh(\omega)\*\mathcal{F}f(x)$. I computed the Fourier transform of $f(x)$, multiplied it with $\sinh(\omega)$ and then computed the inverse Fourier transform.
2. Cauchy integral formula: Generated a matrix $D$ by approximating $\frac{d}{dx}$ using the forward differences and computed the following contour integral: $\sin(D)f=\frac{1}{2\pi i} \oint\limits\_\Gamma \sin(z)(zI-D)^{-1}f \, dz$. Since $D$ is an upper triangular matrix with similar values on the main diagonal, it has one eigenvalue repeated $N$ (size of matrix D) times. Therefore, I chose $\Gamma$ to be a small circle around that eigenvalue.
3. Taylor series: Expanded $\sin(\frac{d}{dx})$ as $\frac{d}{dx}-\frac{d^3}{3!\,dx^3}+\cdots$ and performed discrete differentiations.
4. Scaling and squaring of matrix: Rewrote $\sin(\frac{d}{dx})=\frac{e^{id/dx}-e^{-id/dx}}{2i}$ and used matrix $D$ as defined above.


The calculations in all cases were very unstable. Is there any other (stable) approach that I can try?


Edit: In general, what are the conditions of $g(x)$ and $f(x)$ for $g(\frac{d}{dx})f(x)$ to exist?



How to compute $\sin(\frac{d}{dx})f(x)$?
As noted in the OP, $\sin (d/dx) = (\exp (id/dx) - \exp (-id/dx))/(2i)$, which casts the operator as a combination of two shift operators,
$$
\sin (d/dx) f(x) = \frac{1}{2i} (f(x+i) - f(x-i))
$$
The convergence radius of the Taylor expansion of $f$ around $x$ will have to include $x+i$ and $x-i$.


