# How to implement an implicit iterative method for solving SDEs?

I wish to numerically solve the Black-Scholes SDE as follows $$\begin{array}{lll} dX(t)&=&\mu X(t)dt+\sigma X(t)dW_t, \ \ \ 0\leq t\leq1,\\ X(t_0)&=&X(0), \end{array}$$ with the high-order implicit method of Chaplygin defined by $$\begin{array}{lll} X(0)&=&1., \ \mu=0.9, \ \sigma=2.,\\ X_{n+1}(t)&=&X(0)+\int_0^t{0.9X(s)ds}+\int_0^t{2 X(s)dW_s}\\ &&+\int_0^t{0.9(X_{n+1}(s)-X_n(s))}ds+\int_0^t{2(X_{n+1}(s)-X_n(s))}dW_s. \end{array}$$

This iterative scheme could also be written in its differential form as follows as well

$$\begin{array}{lll} X(0)&=&1,\\ dX_{n+1}(t)&=&0.9X_n(t)dt+2X_n(t)dW_t+0.9(X_{n+1}(t)-X_n(t))dt+2(X_{n+1}(t)-X_n(t))dW_t.\\ \end{array}$$

I have in fact some problem in the implementation of this iterative scheme since it is totally implicit. I will be thankful if anyone could give me some tips/tricks in order to apply it.

Note that Euler-Maruyama method for finding the solution of this SDE has simply been coded in what follows:

ClearAll["Global*"]
SeedRandom[123];
sigma = 2.; mu = 0.9; Xzero = 1.; T = 1.; n = 2^8; dt = T/n;
dW = Sqrt[dt]*RandomVariate[NormalDistribution[0, 1], n];
W = Accumulate[dW];
set = Range[dt, T, dt];
Xtrue = Xzero*Exp[(sigma - 0.5*mu^2)*(set) + mu*W];
NN = Length[set];
solution =
Join[{{0, Xzero}}, Table[{set[[i]], Xtrue[[i]]}, {i, 1, NN}]];
N1 = ListLinePlot[solution, PlotStyle -> {Thin, Blue},
PlotRange -> All];

R = 1; DDt = R*dt; L = n/R; Xtemp = Xzero;
For[ j = 1, j <= L, j++,
{ Winc = Total[Table[dW[[i]], {i, R*(j - 1) + 1, R*j}]];
Xtemp = Xtemp + DDt*sigma*Xtemp + mu*Xtemp*Winc;
X[j] = Xtemp;}];
sol = Flatten[Table[X[j], {j, 1, NN}]];

Numericalsolution =
Join[{{0, Xzero}}, Table[{set[[i]], sol[[i]]}, {i, 1, NN}]];
N2 = ListPlot[Numericalsolution, PlotStyle -> {Thin, Brown, Dashed},
PlotRange -> All];

Show[{N1, N2}, AxesOrigin -> {0, 0}, ImageSize -> 600]
`
-