Previously, the authors derived an analog of the Euler-Maru\-yama method (fEMM) for free stochastic differential equations (fSDEs) and proved strong convergence of order $γ=0.5$ in $L_1(\varphi)$-norm under certain assumptions. In this paper, we study the development of numerical methods for fSDEs which show strong convergence of order $γ=1$ in $L_\infty(\varphi)$. As a side effect, strong convergence of order $γ=0.5$ of fEMM can be extended to $L_p(\varphi)$ for $p\in[1,\infty]$. Utilizing the framework of multiple operator integrals (MOI) we derive a stochastic Itô-Taylor expansion of the solution of the fSDE. It is then possible to identify those free stochastic iterated integrals, which must be discretized in order to obtain strong convergence of order $γ=1$. The non-commutativity imposes additional difficulties showing that the iterated free stochastic integrals can be simulated directly only under special situations, different from the commutative case. We will show, which diffusion terms lead to a Milstein-type method of order $γ=1$. For the cases, where a direct calculation is not possible, we approximate the iterated integrals based on a subdivision of the discretization intervals. As for fEMM, all proposed methods obey strong convergence of order $γ=1$ in $L_p(\varphi),\, 1\leq p\leq \infty$. For all methods developed, we show that the numerical solution is uniformly bounded on finite time intervals.