Friends, this is slightly updated version of the program that corrects for instabilities in the graph that arise from Ito-Taylor tree branches going into negative. Most of the graphs will be perfectly good now. I have also made a few minor changes. Here is the code.

```
Clear[x, t, x0, Y, P1, P2, l1, l2, Z, mu, sigma, mu0, sigma0, alpha, beta, gamma, a0, k, k1, m, L1, L2, YAns, H, Hc, Hz, X1, fX1, DfX1, Nn, n, HCoeff, ZCoeff, DZCoeff, zz, KIndex, MIndex, MaxIndex, Y1, Y2, CoeffDX1, CoeffDX2, CoeffdtIntegral, CoeffdzIntegral, qIndex1, qIndex21, qIndex22, qIndex23, qIndex0, qqIndex1, qqIndex2, n1, m1, l3, P, k0, m0];
P1 = 4;
P2 = 4;
P3 = 4;
L1 = 0;
L2 = 0;
P = If[P1 > P2, P1, P2];
Array[Y, {P + L1 + 1, P + L1 + 1, P + L2 + 1}, {0, 0, 0}];
For[k = 0, k <= P + L1, k++, (For[m = 0, m <= P + L1, m++, (For[n = 0, n <= P + L2, n++, (Y[k, m, n] = 0);]);]);];
Array[KIndex, {P + 1}, {0}];
For[k = 0, k <= P + 1, k++, If[k > 0, (k0 = k - 1;
KIndex[k] = (k0*(k0 + 1)*(2*k0 + 1) + 9*k0*(k0 + 1) + 12*k0)/12 + 1), KIndex[k] = 0]];
Array[MIndex, {P1 + 1}, {0}];
For[m = 0, m <= P + 1, m++, If[m > 0, (m0 = m - 1; MIndex[m] = m0*(m0 + 1)/2 + m0 + 1),
MIndex[m] = 0]];
MaxIndex = 1 + 2*MIndex[Round[P/2 + 1]]*KIndex[Round[P/2 + 1]];
Array[Y1, MaxIndex, 0];
Array[Y2, MaxIndex, 0];
For[m = 0, m < MaxIndex, m++, (Y1[m] = 1; Y2[m] = 0;)];
Array[CoeffdtIntegral, {P + 1, P + 1}, {0, 0}];
Array[CoeffdzIntegral, {P + 1, P + 1}, {0, 0}];
For[k = 0, k < P, k++, (For[m = 0, m <= k, m++, (CoeffdtIntegral[k - m, m] = 1/(m + 1)*(1 - Sqrt[k - m]/Sqrt[2*m + (k - m) + 2]);
CoeffdzIntegral[k - m, m] = 1/Sqrt[(2*m + (k - m) + 1)*(k - m + 1)];)])];
a0 = 1;
mu0 = .2;
sigma0 = 1.0;
alpha = 1;
gamma = .65;
beta = 1;
Y[L1, L1, L2] = a0;
For[k = 0, k < P, k++, (For[m = 0, m <= k, m++, (l3 = L2 + k - m;
For[n = 0, n <= m, n++, (l1 = L1 + m - n; l2 = L1 + n; qIndex1 =
MIndex[m] + n + 1; QIndex1 = (qIndex1 - 1)*KIndex[P - k + 1];
qIndex21 = MIndex[m + 1] + n + 1;
QIndex21 = (qIndex21 - 1)*KIndex[P - k];
qIndex22 = MIndex[m + 1] + (n + 1) + 1;
QIndex22 = (qIndex22 - 1)*KIndex[P - k];
qIndex23 = MIndex[m] + n + 1;
QIndex23 = (qIndex23 - 1)*KIndex[P - k];
CoeffDX1 = (alpha + (l1 - L1)*beta + (l2 - L1)*2*gamma + (l3 - L2)*gamma - (l1 - L1) - 2*(l2 - L1) - (l3 - L2));
CoeffDX2 = CoeffDX1 - 1; MuDX1 = CoeffDX1*mu0;
SigmaDX1 = CoeffDX1*sigma0;
Sigma2DX2 = .5*CoeffDX1*CoeffDX2*sigma0^2;
For[k1 = 0, k1 < P - k, k1++, (For[m1 = 0, m1 <= k1, m1++, (For[n1 = 0, n1 <= m1, n1++, (qqIndex2 = KIndex[k1] + MIndex[m1] + n1 + 1;
qqIndex1 = KIndex[k1 + 1] + MIndex[m1 + 1] + n1 + 1;
Y2[qqIndex2 + QIndex21] = Y2[qqIndex2 + QIndex21] + Y1[qqIndex1 + QIndex1]*CoeffdtIntegral[k1 - m1, m1]*MuDX1; qqIndex1 =
KIndex[k1 + 1] + MIndex[m1 + 1] + (n1 + 1) + 1;
Y2[qqIndex2 + QIndex22] = Y2[qqIndex2 + QIndex22] + Y1[qqIndex1 + QIndex1]*CoeffdtIntegral[k1 - m1, m1]*Sigma2DX2;
qqIndex1 = KIndex[k1 + 1] + MIndex[m1] + n1 + 1;
Y2[qqIndex2 + QIndex23] = Y2[qqIndex2 + QIndex23] + Y1[qqIndex1 + QIndex1]*CoeffdzIntegral[k1 - m1, m1]*SigmaDX1;)];)])];
Y[l1 + 1, l2, l3] = Y2[1 + QIndex21];
Y[l1, l2 + 1, l3] = Y2[1 + QIndex22];
Y[l1, l2, l3 + 1] = Y2[1 + QIndex23];)];)];
For[qIndex0 = 0, qIndex0 < MaxIndex, qIndex0++, (Y1[qIndex0] = Y2[qIndex0]; Y2[qIndex0] = 0;)];)];
Nn = 81;
P2 = 4;
Array[HCoeff, {P2 + L2 + 1, Nn}, {0, 0}];
Array[Hz, {P2 + L2 + 1, P2 + L2 + 1}, {0, 0}];
For[k = 0, k < P2 + L2 + 1, k++,
For[m = 0, m < P2 + L2 + 1, m++, Hz[k, m] = 0;]];
Hz[0, 0] = 1.0;
For[k = 0, k < P2 + L2, k++, (For[m = 0, m <= k, m++, (Hz[k + 1, m + 1] = Hz[k, m];
If[m >= 1, Hz[k + 1, m - 1] = Hz[k + 1, m - 1] - m*Hz[k, m], 0];)])];
x0 = 1;
t = .1250;
Tt = 8;
For[k = 0, k < P2 + L2 + 1, k++, (For [nn = 0, nn < Nn, nn++, HCoeff[k, nn] = 0;])];
Array[x, Nn, 0];
For[nn = 0, nn < Nn, nn++, x[nn] = x0];
Array[ZCoeff, {P2 + L2 + 1, Nn}, 0];
Array[X1, Nn, 0];
Array[DfX1, Nn, 0];
Array[fX1, Nn, 0];
Array[zz, Nn, 0];
For[nn = 0, nn < Nn, nn++, (X12[nn] = 0; DfX12[nn] = 0;)];
Array[DZCoeff, {P2 + L2, Nn}, {0, 0}];
For[tt = 0, tt < Tt, tt++, (For[nn = 0, nn < Nn, nn++, (For[k = 0, k <= P3, k++, HCoeff[k, nn] = 0];
For[k = 0, k <= P3, k++, (For[m = 0, m <= If[k < P1, k, P1], m++, (l0 = L1 + m; l3 = L2 + k - m;
For[n = 0, n <= If[m < P2, m, P2], n++, (l1 = L1 + m - n; l2 = L1 + n;
HCoeff[l3, nn] = HCoeff[l3, nn] + Y[l1, l2, l3]*x[nn]^(alpha + (l1 - L1)*beta + (l2 - L1)*2*gamma + (l3 - L2)*gamma - (l1 - L1) - 2*(l2 - L1) - (l3 - L2))*t^(l1 + l2 - L1 + .5*l3);)];)])])];
For[k = 0, k < P2 + L2 + 1, k++, For[nn = 0, nn < Nn, nn++, ZCoeff[k, nn] = 0]];
For[nn = 0, nn < Nn, nn++, (For[k = 0, k < P2 + L2 + 1, k++,
For[m = 0, m <= k, m++, ZCoeff[m, nn] = ZCoeff[m, nn] + Hz[k, m]*HCoeff[k, nn]]];)];
For[nn = 0, nn < Nn, nn++, (X12[nn] = ZCoeff[0, nn]; zz[nn] = nn*.1 - 4;
For[k = 1, k < P2 + L2 + 1, k++, X12[nn] = X12[nn] + ZCoeff[k, nn]*(1/(Tt)^(.5)*zz[nn])^k;]; x[nn] = X12[nn];)];)];
CFlag = 1;
For[nn = Nn - 1, nn > 1, nn--, (If[((X12[nn] > X12[nn - 1]) && (CFlag == 1)), (X13[nn] = X12[nn]; nnStart = nn;), (X13[nn] = 0; CFlag = 0;)])];
nnStart
For[nn = 1, nn < Nn - 1, nn++, zz[nn] = (nn*.1 - 4);];
For[nn = nnStart, nn < Nn - 1, nn++, (
DfX12[nn] = (X12[nn + 1] - X12[nn - 1])/(zz[n + 1] - zz[n - 1]);)];
For[nn = nnStart, nn < Nn, nn++, (zz[nn] = nn*.1 - 4;
fX12[nn] = PDF[NormalDistribution[0, 1], zz[nn]]/Abs[DfX12[nn]];)];
data12 = Table[{X12[nn], fX12[nn]}, {nn, nnStart, 60}];
p12 = ListLinePlot[data12]
```

Here is the output graph of the Ito-Taylor tree program overlaid with graph from true solution using analytic formula.