  In this paper, we develop an accurate pseudospectral method to approximate
numerically the Riesz-Feller operator $D_\gamma^\alpha$ on $\mathbb R$, where
$\alpha\in(0,2)$, and $|\gamma|\le\min\{\alpha, 2 - \alpha\}$. This operator
can be written as a linear combination of the Weyl-Marchaud derivatives
$\mathcal{D}^{\alpha}$ and $\overline{\mathcal{D}^\alpha}$, when
$\alpha\in(0,1)$, and of $\partial_x\mathcal{D}^{\alpha-1}$ and
$\partial_x\overline{\mathcal{D}^{\alpha-1}}$, when $\alpha\in(1,2)$.
  Given the so-called Higgins functions $\lambda_k(x) = ((ix-1)/(ix+1))^k$,
where $k\in\mathbb Z$, we compute explicitly, using complex variable
techniques, $\mathcal{D}^{\alpha}[\lambda_k](x)$,
$\overline{\mathcal{D}^\alpha}[\lambda_k](x)$,
$\partial_x\mathcal{D}^{\alpha-1}[\lambda_k](x)$,
$\partial_x\overline{\mathcal{D}^{\alpha-1}}[\lambda_k](x)$ and
$D_\gamma^\alpha[\lambda_k](x)$, in terms of the Gaussian hypergeometric
function ${}_2F_1$, and relate these results to previous ones for the
fractional Laplacian. This enables us to approximate
$\mathcal{D}^{\alpha}[u](x)$, $\overline{\mathcal{D}^\alpha}[u](x)$,
$\partial_x\mathcal{D}^{\alpha-1}[u](x)$,
$\partial_x\overline{\mathcal{D}^{\alpha-1}}[u](x)$ and
$D_\gamma^\alpha[u](x)$, for bounded continuous functions $u(x)$. Finally, we
simulate a nonlinear Riesz-Feller fractional diffusion equation, characterized
by having front propagating solutions whose speed grows exponentially in time.
