I am posting a new version of the old programs correcting the error. I have run the new program on density of square of a lognormal SDE with drift. And it can then be compared to analytical solution for the density of squared lognormal. Here is the new code with Ito-Taylor.

`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];`

P1=24;

P2=24;

L1=0;

L2=0;

Array[Y,{P1+L1+1,P1+L1+1,P2+L2+1},{0,0,0}];

For[k=0,k<=P1+L1,k++,(For[m=0,m<=P1+L1,m++,(For[n=0,n<=P2+L2,n++,(Y[k,m,n]=0);]);]);];

a0=1;

mu0=.25;

sigma0=.75;

alpha=2;

gamma=1.0;

beta=1;

Y[L1,L1,L2]=a0;

For[k=0,k<P1,k++,(For[m=0,m<=If[k<=P2,k,P2],m++,(l0=L1+k-m;l3=L2+m;

For[n=0,n<=l0-L1,n++,(l1=l0-n;l2=L1+n;

Y[l1+1,l2,l3]=Y[l1+1,l2,l3]+SetPrecision[(1.0/(l1+l2+1)*(1.0-Sqrt[l3]/Sqrt[2*l1+2*l2+l3+2]))*Y[l1,l2,l3]*(alpha+(l1-L1)*beta+(l2-L1)*2*gamma+(l3-L2)*gamma-(l1-L1)-2*(l2-L1)-(l3-L2))*mu0,50];

Y[l1,l2+1,l3]=Y[l1,l2+1,l3]+SetPrecision[(1.0/(l1+l2+1)*(1.0-Sqrt[l3]/Sqrt[2*l1+2*l2+l3+2]))*Y[l1,l2,l3]*.5*(alpha+(l1-L1)*beta+(l2-L1)*2*gamma+(l3-L2)*gamma-(l1-L1)-2*(l2-L1)-(l3-L2))*(alpha+(l1-L1)*beta+(l2-L1)*2*gamma+(l3-L2)*gamma-(l1-L1)-2*(l2-L1)-(l3-L2)-1)*sigma0^2,50];

If[l3<L2+P2,(Y[l1,l2,l3+1]=Y[l1,l2,l3+1]+SetPrecision[(1.0/Sqrt[2*l1+2*l2+l3+1]/Sqrt[l3+1])*Y[l1,l2,l3]*(alpha+(l1-L1)*beta+(l2-L1)*2*gamma+(l3-L2)*gamma-(l1-L1)-2*(l2-L1)-(l3-L2))*sigma0,50]),0];);];)];);];

Array[HCoeff,P2+L2+1,0];

For[k=0,k<P2+L2+1,k++,HCoeff[k]=0;]

x=1;

t=2;

For[k=0,k<=P1,k++,(For[m=0,m<=If[k<=P2,k,P2],m++,(l0=L1+k-m;l3=L2+m;

For[n=0,n<=l0-L1,n++,(l1=l0-n;l2=L1+n;

HCoeff[l3]=HCoeff[l3]+SetPrecision[Y[l1,l2,l3]*x^(alpha+(l1-L1)*beta+(l2-L1)*2*gamma+(l3-L2)*gamma-(l1-L1)-2*(l2-L1)-(l3-L2))*(Sqrt[t^(2*l1+2*l2-2*L1+l3)]),50];)];)];);];

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]=SetPrecision[Hz[k,m],50];

If[m>=1,Hz[k+1,m-1]=SetPrecision[Hz[k+1,m-1]-m*Hz[k,m],50],0];)])];

Array[ZCoeff,P2+L2+1,0];

For[k=0,k<P2+L2+1,k++,ZCoeff[k]=0];

For[k=0,k<P2+L2+1,k++,(For[m=0,m<=k,m++,ZCoeff[m]=ZCoeff[m]+SetPrecision[Hz[k,m]*HCoeff[k],50]])];

Array[DZCoeff,P2+L2,0];

For[k=0,k<P2+L2,k++,DZCoeff[k]=0];

For[k=0,k<P2+L2,k++,DZCoeff[k]=SetPrecision[ZCoeff[k+1]*(k+1),50]];

Nn=81;

Array[X1,Nn,0];

Array[DfX1,Nn,0];

Array[fX1,Nn,0];

Array[zz,Nn,0];

For[n=0,n<Nn,n++,(X1[n]=0;DfX1[n]=0;fX1[n]=0;zz[n]=0;)];

For[n=0,n<Nn,n++,(zz[n]=n*.1-4;X1[n]=0;DfX1[n]=0;fX1[n]=0;

X1[n]=ZCoeff[0];DfX1[n]=DZCoeff[0];For[k=1,k<P2+L2+1,k++,(X1[n]=X1[n]+SetPrecision[ZCoeff[k]*zz[n]^k,50];If[k<P2+L2,DfX1[n]=DfX1[n]+SetPrecision[DZCoeff[k]*zz[n]^k,50],0];)];

fX1[n]=SetPrecision[PDF[NormalDistribution[0,1],zz[n]]/Abs[DfX1[n]],50];)];

data1=Table[{X1[n],fX1[n]},{n,0,58}];

p1=ListLinePlot[data1]

And here is the code for analytical density of square of lognormal SDE with drift.

I hope that this new program will be helpful.

Here is the graph comparison image. The parameters are mu=.25, sigma=.75, X0=1,t=2 for geometric Brownian motion/lognormal SDE. And graphs show the density of squared X(t) that is X(t)^2. The following graph shows the overlaid comparison for analytical solution with Ito-Taylor expansion. Both images are indistinguishable. I have removed a lot of the area from the tail where probability was miniscule to focus on the main area.

Here is graph comparison for the lognormal SDE density with drift but that is not squared.

It is again overlaid graph comparison with the same parameters as described earlier but it is the SDE variable and not its square. Here is the overlaid comparison graph. I have again tried to focus on the main graph body and not on the tails.

Here is another overlaid graph comparison of ito-taylor expansion and analytical formula for five year density of lognormal SDE with slightly decreased volatility. Here mu=.25, sigma=.65,t=5, X0=1.