# Improving the speed on an iterated differential system

I'm working on a simulation where I need to solve a differential system for a function, and then use the solution to this system to calculate another function that is used as an input to this differential system. This has to be iterated several times.

In code, the situation looks something like this:

solvefunction[fun_] :=
NDSolveValue[{g'[x] == fun[x]*g[x], g == 1}, {g}, {x, 0, 1}]


This calculates the first function as a solution to a system (here simplified to only one equation). Then, the solution to this system is used to build a new function:

findfunction[fun_] := solvefunction[fun[1 - #^2] &][]


Which then has to be used as the input to the previous system. I can iterate using:

Nest[findfunction, #&, 4]


The first iterations are fast but for n>4 it becomes very slow. I suspect this recurrence is not very efficient from a memory point of view. Is there a generic way to achieve this result in a faster way ?

• I have not tried it but perhaps ParametricNDSolve is useful. – user21 Sep 7 '16 at 15:25
• It looks like the 4th one blows up real good (i.e. it's equal to 10^10 at x = 1), and so perhaps the values of the function in the next iteration get too large for NDSolve to handle. Can you scale the functions back to having a value of 1 at x = 1 at each step? It seems to me that that won't change the functional form of each iteration. – march Sep 7 '16 at 15:49
• Hmm. Maybe I'm wrong about it not changing the form. I'm not sure how to fix that. But after experimenting, I'm pretty sure the problem is that the function values are blowing up, as I described above. – march Sep 7 '16 at 15:59
• You are right, that seems to be a problem. By changing the differential equation to g'[x] == -fun[x]*g[x], it does compute a lot faster. – Whelp Sep 8 '16 at 9:55