Simulating a Stochastic Integral of OU process The stochastic integral I want to simulate is 
$$\int_{0}^{1}J_c(s)dJ_c(s)$$
where $J_c(s) = \int_{0}^{s}e^{-c(s-r)}dB(r)$, is an OU process.

I simulate the data using Matlab and the sample codes are as follows. n is the sample size.

dB_x = mvnrnd(zeros(K, 1), Omega_xx/n, n); % n-by-m matrix
B_x = cumsum(dB_x,1); % n-by-m matrix
B_x = [zeros(1,K); B_x];
J_c = exp(-C*[1:n]/n)'.*cumsum(exp(C*[1:n]/n).*[zeros(K,1) dB_x(1:n-1,:)'])';
dJ_c=[J_c(1,:); diff(J_c)]; % m-by-n matrix
int_J =  1/n*J_c'*dJ_c;


I am not sure whether my calculation is correct or not, can anybody give me a correct way of simulating this integral?
Note that
\begin{align*}
dJ_c(s) = dB_s -cJ_c(s) ds.
\end{align*}
Moreover, $d\langle J_c, J_c\rangle_t = dt$.
Then, by Ito's lemma,
\begin{align*}
\int_0^1J_c(s) dJ_c(s) &amp;= \int_0^1 d\left(J_c^2(s) -s \right)\\
&amp;= J_c^2(1) - 1.
\end{align*}
Since
\begin{align*}
J_c(1) &amp;= \int_0^1 e^{-c(1-s)}dB_s
\end{align*}
is normal with zero mean and the variance given by
\begin{align*}
E\left(J_c^2(1) \right) &amp;= \int_0^1e^{-2c(1-s)} ds\\
&amp;=\frac{1}{2c}(1-e^{-2c}).
\end{align*}
Then \begin{align*}
J_c(1) = \sqrt{\frac{1}{2c}(1-e^{-2c})}\, Z,
\end{align*}
and
\begin{align*}
\int_0^1J_c(s) dJ_c(s) &amp;= \frac{1}{2c}(1-e^{-2c})\, Z^2 -1.
\end{align*}
where $Z$ is a standard normal random variable. The simulation becomes a draw of a sample from a standard normal random variable.
Finally, I found it wrong. The correct way is

dB_x = mvnrnd(zeros(K, 1), Omega_xx/n, n); % n-by-m matrix
B_x = cumsum(dB_x,1); % n-by-m matrix
B_x = [zeros(1,K); B_x];
J_c = exp(C*[1:n]/n)'.*cumsum(exp(-C*[1:n]/n).*[zeros(K,1) dB_x(1:n-1,:)'])';
dJ_c=[J_c(1,:); diff(J_c)]; % m-by-n matrix
int_J = J_c'*dJ_c;