  We examine the numerical approximation of a quasilinear stochastic
differential equation (SDE) with multiplicative fractional Brownian motion. The
stochastic integral is interpreted in the Wick-It\^o-Skorohod (WIS) sense that
is well defined and centered for all $H\in(0,1)$. We give an introduction to
the theory of WIS integration before we examine existence and uniqueness of a
solution to the SDE. We then introduce our numerical method which is based on
the theoretical results in \cite{Mishura2008article, Mishura2008} for $H\geq
\frac{1}{2}$. We construct explicitly a translation operator required for the
practical implementation of the method and are not aware of any other
implementation of a numerical method for the WIS SDE. We then prove a strong
convergence result that gives, in the non-autonomous case, an error of
$O(\Delta t^H)$ and in the non-autonomous case $O(\Delta t^{\min(H,\zeta)})$,
where $\zeta$ is a H\"older continuity parameter. We present some numerical
experiments and conjecture that the theoretical results may not be optimal
since we observe numerically a rate of $\min(H+\frac{1}{2},1)$ in the
autonomous case. This work opens up the possibility to efficiently simulate
SDEs for all $H$ values, including small values of $H$ when the stochastic
integral is interpreted in the WIS sense.
