OK, a source for SDEs has been found; have you a similar source for ODEs?

Iterated integrals can be used in so many different ways and that does not mean that everywhere iterated integrals are used, it is the same setting in which I used them for ODEs. The PDF you have posted has little relevance to my work on ODEs. Can you please suggest something where iterated integrals have been used explicitly for ODEs in a similar manner to how I used them? 
I developed the idea explicitly in the case of 1st order, 2nd order and nth order ODEs case and also presented programs to solve systems of  ODEs of nth order.

OK. The integral and ODE approach are mirror images of each other.
It is really up to you to show it is original work. A good idea is to locate ODE academics and ask them.
Iterated integrals is not my cup of tea, since I prefer numerics. Each to his own, I suppose.

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.

Code: Select all

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]

I will take a few more days to complete the paper but here are a few formulas for friends. Here I give the general algorithm for the type of integrals where only dz and dt integrals follow each other in any order. Like many other possibilities, one four order repeated integral with two dz uncertainties could be[$]\int _0^t \text{$\sigma $dz}(s)\int _0^s\text{du}\int _0^u\text{dv}\int _0^v\sigma  \text{dz}(w)[$]The value of the above kind of integrals, as long as there is no explicit t, does not depend upon the sequence order(whether dz integral comes first or dt integral comes first does not change the value of the total integrals) so variance can be easily calculated by Ito Isometry. The order of hermite polynomial depends upon the number of normal uncertainties so is equal to number of dz integrals. Here is the simple formula to calculate the hermite polynomial representation of the integrals which equals[$] H_n(N)*\frac{\sqrt{\text{Variance}}}{\sqrt{n!}}[$]here n, which is the order of hermite polynomials, equals the number of normal uncertainties.   I will give a few examples[$]\int _0^t \text{$\sigma $dz}(s)\int _0^s\text{du}=H_1(N)\frac{\text{$\sigma $t}^{1.5}}{\sqrt{1!}\sqrt{3}}[$][$]\int _0^t \text{$\sigma $dz}(s)\int _0^s\text{$\sigma $dz}(u)=H_2(N)\frac{\sigma ^2t}{\sqrt{2!}\sqrt{2}}[$][$]\int _0^t \text{$\sigma $dz}(s)\int _0^s\text{du}\int _0^u\text{$\sigma $dz}(v)=\int _0^t \text{$\sigma $dz}(s)\int _0^s\text{$\sigma $dz}(u)\int _0^u\text{dv}=\int _0^t \text{ds}\int _0^s\text{$\sigma $dz}(u)\int _0^u\text{$\sigma $dz}(v)=H_2(N)\frac{\sigma ^2t^2}{\sqrt{2!}\sqrt{12}}[$][$]\int _0^t \text{$\sigma $dz}(s)\int _0^s\text{du}\int _0^u\text{dv}=\int _0^t \text{ds}\int _0^s\text{$\sigma $dz}(u)\int _0^u\text{dv}=\int _0^t \text{ds}\int _0^s\text{du}\int _0^u\text{$\sigma $dz}(v)=H_1(N)\frac{\sigma ^2t^{2.5}}{\sqrt{1!}\sqrt{20}}[$][$]\int _0^t \text{$\sigma $dz}(s)\int _0^s\text{$\sigma $dz}(u)\int _0^u\text{$\sigma $dz}(v)=H_3(N)\frac{\sigma ^3t^{1.5}}{\sqrt{3!}\sqrt{6}}[$]Fourth order and higher integrals can be evaluated to perfect precision using this method so I am not giving examples for that.