Does anybody know how to calculate self-calling differential operators fast?

Hello friends,

I have a code with a for-loop where inside the for-loop there is a differential operator which uses a function, say A0, as input and performs

some operations and fthe results of the calculations are asigned to A0 again to be used for the next iteration. Intuitively, we understand that

the expressions in the for-loop get more and more complex at each iteration so that the computation time increases.

Long story short: In bellow I attach my code. The double for-loop (with tic and toc) becoms computationlly expensive

very when I increase J and K. I hope you have an idea to imrpve the speed of this code!

Thanks lot in advance!

Babak

clc;
J=4;K=4;
H=sym('H',[J 1]);dt=0.1;
syms x par1 par2 par3 par4 par5
for j=1:J
    H(j)=horner(simplify(exp(x^2/2)*diff(exp(-x^2/2),j)),x);
end
syms y y0 positive
muY(y)=-(y*(y*(y*(par1*y*par5^2 - par1*par2*par5) + par2*par3 + par1*par5^2) - par1*par2*par5))/(par2*par5*y^2 + par2*par5);
Hz=subs(H,x,(y-y0)/dt^(1/2));
eta=sym('eta',[J 1]);
tic;
for j=1:J
    A0=Hz(j);S=A0;
    for k=1:K
        A0=muY(y)*diff(A0,y)+1/2*diff(A0,y,2);
        S=S+A0*dt^k/factorial(k);
    end
    eta(j)=1/factorial(j)*limit(S,y,y0);    
end 
toc; 

[Walter Roberson]

I wonder...

If you were to use

syms MUY(y)

and

A0=MUY(y)*diff(A0,y)+1/2*diff(A0,y,2);

then you should be able to get out a response in terms of powers of MUY(y) and derivatives of MUY(y) . At that point, you can use findSymType() to locate derivatives of MUY, and find the maximum derivative, as well as building up a list of which derivatives are used. Loop through iteratively building up the derivatives to the maximum used. Then subs() for the derivative actually used in S.

I think doing this should lower the cost of doing the derivatives of muY. The subs() at the end might be expensive, but I suspect you will still end up saving time.