| T=1; NN=10; h=T/NN; R=10000; q=0; |
| for r=1:R, x=2; y=0; xx=2; yy=0; |
| for m=1:NN, hh=h/2; N=10; |
| [wL,vL,J1L,J2L] = F_year_J(N,hh); |
| = F_year_J(N,hh); |
| w=wL+wr; v=vL+vr; J1=J1L+J1r+wrvL; J2=J2L+J2r+wLvr; |
| u=x+ywL+(x+(m−1)h)vL+exp(−y. 2)(1/2)(wL. 2−hh) |
| +(x−y)J1L+yJ2L+(x+(m−1)h)(1/2)(vL. 2−hh); |
| y=y+exp(−y. 2)wL+(x−y)vL−2yexp(−2y. 2)(1/2)(wL. 2−hh) |
| +(y−exp(−y. 2))J2L+(y+(m−1)h)(1/2)(vL. 2−hh) |
| −2yexp(−y. 2)(x−y)J1L; x=u; |
| u=x+ywr+(x+(m−1/2)h)vr+exp(−y. 2)(1/2)(wr. 2−hh) |
| +(x−y)J1r+yJ2r+(x+(m−1/2)h)(1/2)(vr. 2−hh); |
| y=y+exp(−y. 2)wr+(x−y)vr−2yexp(−2y. 2)(1/2)(wr. 2−hh) |
| +(y−exp(−y. 2))J2r+(y+(m−1/2)h)(1/2)(vr. 2−hh) |
| −2yexp(−y. 2)(x−y)J1r; x=u; |
| u=xx+yyw+(xx+(m−1)h)v+exp(−yy. 2)(1/2)(w. 2−h) |
| +(xx−yy)J1+yyJ2+(xx+(m−1)h)(1/2)(v. 2−h); |
| yy=yy+exp(−yy. 2)w+(xx−yy)v−2yyexp(−2yy. 2)(1/2)(w. 2−h) |
| +(yy−exp(−yy. 2))J2+(yy+(m−1)h)(1/2)(v. 2−h) |
| −2yyexp(−yy. 2)(xx−yy)J1; xx=u; end |
| q=q+abs(x−xx)+abs(y−yy); end |
| (q/R) |
