Here is a program that will calculate the density of an SDE using hermite polynomials. In the following program, I have used 24 hermite polynomials for expansion of driftless CEV SDE density at two years. I have kept the drift at zero and this is just an initial experimentation. I have not carefully checked everything and posting the program in the hope that it would be helpful for others. I will stay away from this project for a few days and possibly come back next week and make more changes. Here is the code.

`In[83]:= 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,Hz,X1,fX1,DfX1,Nn,n,HCoeff,ZCoeff,DZCoeff,zz];`

In[84]:= P1=24;

In[85]:= P2=24;

In[86]:= L1=0;

In[87]:= L2=0;

In[88]:= Array[Y,{P1+L1+1,P1+L1+1,P2+L2+1},{0,0,0}];

In[89]:= 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);]);]);];

In[90]:= a0=1;

In[91]:= mu0=0;

In[92]:= sigma0=.75;

In[93]:= alpha=1;

In[94]:= gamma=.85;

In[95]:= beta=0;

In[96]:= Y[L1,L1,L2]=a0;

In[97]:= 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[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[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[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];);];)];);];

In[98]:= Array[HCoeff,P2+L2+1,0];

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

In[100]:= x=1;

In[101]:= t=2;

In[102]:= 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)]*Sqrt[(2*l1+2*l2-2*L1)!])/(Sqrt[(2*l1+2*l2-2*L1+l3)!] *((l1+l2-L1)!))/Sqrt[l3!],50];)];)];);];

In[103]:= Array[Hz,{P2+L2+1,P2+L2+1},{0,0}];

In[104]:= For[k=0,k<P2+L2+1,k++,For[m=0,m<P2+L2+1,m++,Hz[k,m]=0;]];

In[105]:= Hz[0,0]=1.0;

In[106]:= 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];)])];

In[107]:= Array[ZCoeff,P2+L2+1,0];

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

In[109]:= 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]])];

In[110]:= Array[DZCoeff,P2+L2,0];

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

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

In[113]:= Nn=81;

In[114]:= Array[X1,Nn,0];

In[115]:= Array[DfX1,Nn,0];

Array[fX1,Nn,0];

Array[zz,Nn,0];

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

In[119]:= 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];

)];

In[120]:= data1=Table[{X1[n],fX1[n]},{n,0,80}];

In[122]:= ListLinePlot[data1]