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.