We approximate the solution of some linear systems of SDEs driven by a fractional Brownian motion $B^H$ with Hurst parameter $H\in(\frac{1}{2},1)$ in the Wick--Itô sense, including a geometric fractional Brownian motion. To this end, we apply a Donsker-type approximation of the fractional Brownian motion by disturbed binary random walks due to Sottinen. Moreover, we replace the rather complicated Wick products by their discrete counterpart, acting on the binary variables, in the corresponding systems of Wick difference equations. As the solutions of the SDEs admit series representations in terms of Wick powers, a key to the proof of our Euler scheme is an approximation of the Hermite recursion formula for the Wick powers of $B^H$.