In the series expansion to obtain the form of higher order Milstein schemes you need to evaluate these two sister integrals (where W_t is the Wiener process): $$ \int_{t_i}^{t_{i+1}}dW_u\int_{t_i}^{u}ds\\ \int_{t_i}^{t_{i+1}}du\int_{t_i}^{u}dW_s, $$ after I evaluate the innermost integral I obtain respectively: $$ \int_{t_i}^{t_{i+1}}dW_u(u-t_i)\\ \int_{t_i}^{t_{i+1}}du (W_u-W_{t_i}). $$ However I find myself stuck at this point since I cannot explicitly integrate the Wiener process under $du$ and vice-versa $u$ under $dW_u$. How can I solve these two integrals?
P.S. I was able to solve for the sum of the two integrals $$ \int_{t_i}^{t_{i+1}}dW_u\int_{t_i}^{t_{i+1}}ds+\int_{u}^{t_{i+1}}du\int_{t_i}^{u}dW_s=(t_{i+1}-t_i)^{3/2}\zeta; $$ where $\zeta$ is a Gaussian random variable however I could not use this result in any way.
P.P.S Do you have any good references for these kind of integrals? I am basing my question on the second chapter of the Kloeden&Platten book but they provide little insight on the calculation.
In a numerical solver all you need to know are the variances of those sister integrals and their covariance. From the Ito isometry we get \begin{align}\tag{1} \mathbb E\Big[\Big(\textstyle\int_a^bu-a\,dW_u\Big)^2\Big]=\int_a^b(u-a)^2\,du=\frac{(b-a)^3}{3}\,. \end{align} Next, using integration-by-parts, \begin{align}\tag{2} \textstyle\int_a^b W_udu&=-\textstyle\int_a^bu\,dW_u+b\,W_b-a\,W_a \end{align} so that \begin{align}\require{cancel} \textstyle\int_a^b W_u-W_a\,du&=-\textstyle\int_a^bu\,dW_u+b\,W_b-\cancel{a\,W_a}-b\,W_a+\cancel{a\,W_a}\\[2mm] &=\textstyle\int_a^bb-u\,dW_u\,.\tag{3} \end{align} Similarly to (1) the variance of this is $$\tag{4} \mathbb E\Big[\Big(\textstyle\int_a^bb-u\,dW_u\Big)^2\Big]=\int_a^b(b-u)^2\,du=\frac{(b-a)^3}{3}\,. $$
The covariance is \begin{align}\tag{5} \mathbb E\Big[\Big(\textstyle\int_a^bu-a\,dW_u\Big)\Big(\textstyle\int_a^bb-u\,dW_u\Big)\Big]=\int_a^b(u-a)(b-u)\,du=\textstyle\frac{(b-a)^3}{6}\,. \end{align}
The standard reference is Kloeden & Platen, Numerical Solution of Stochastic Differential Equations.