# Using output of NDSolve in further NDSolve calculations [closed]

I am trying to use the interpolating functions found using NDSolve to solve a further differential equation where the first solutions depend on two arbitrary parameters.

As a simple example, suppose I have two coupled ODES that I can write as:

func[r_,u_]:=y''[x] == u*y'[x]+r^2*y[x]*x^3
sol[r_,u_]:= NDSolve[{func[r,u],y[0]==r,y'[0]==r*u},y,{x,0,50}]
func2[r_,u_]:= (z''[x]==y'[x]*r+u*y[x]*z'[x])/.sol[r,u]
sol2[r_,u_]:= NDSolve[{func2[r,u],z[0]==1,z'[0]==0},z,{x,0,50}]


The problem is that I get the error message when running the fourth line telling me that replacing y[x] etc. is not a good replacement rule.

In reality, my system is a horrible non-linear set f coupled ODEs with one final equation that depends on the solutions of the others. Furthermore, r and u also control the boundary conditions but I think this example suffices because I think the problems are caused by the fact that there is r and u dependence in the ODE.

I can see many posts where the output of NDSolve does not depend on any parameters, in which case it seems trivial to just define a function that does what I want e.g.

f1[x_]:=Evaluate[First[y[x]/.NDSolve[{y[x]''+y[x]'+x^3==0,y[0]==1,y'[0]==0},y,{x,0,50}]]]

sol=NDSolve[{z'[x]==^3*y[x]+y'[x],z[0]=0},z,{x,0,50}]


would work if it were not for u.

I am also aware that I can solve both equations simultaneously but the second equation is stiff so I would like to use a different method for it. Actually, the first set are DAE's whilst the last one is not so splitting the two would help me a lot.

## closed as unclear what you're asking by Michael E2, m_goldberg, MarcoB, ubpdqn, Yves KlettAug 3 '16 at 5:47

Please clarify your specific problem or add additional details to highlight exactly what you need. As it's currently written, it’s hard to tell exactly what you're asking. See the How to Ask page for help clarifying this question. If this question can be reworded to fit the rules in the help center, please edit the question.

• How are you calling the function? – Feyre Aug 1 '16 at 15:29
• Welcome! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour and check the faqs! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! – user9660 Aug 1 '16 at 15:33
• – user9660 Aug 1 '16 at 15:34
• Thanks for the help guys. I actually wrote this code in by hand because my actual code is a lot more complicated and my real ODE is very long. Feyre, what do you mean how and I calling the function? – jsaxon Aug 1 '16 at 15:41
• sol2[r_,u_]:=  is set delayed, so I assume you are calling it in some form of sol2[1,2] – Feyre Aug 1 '16 at 16:09

Well, I made some changes, but then realized I hadn't changed anything important. I shortened the integration to {x, 0, 2} to save time. The original code works with this change. I memoized sol so that it's computed only once for each {r, u}. You stick First@ on NDSolve[] to remove a set of braces, but it's not important.

Clear[func, func2, sol, sol2, x, y, r, u];

func[r_, u_] := y''[x] == u*y'[x] + r^2*y[x]*x^3;
sol[r_?NumericQ, u_?NumericQ] :=
sol[r, u] = (* saves the value so that it is computed once *)
NDSolve[{func[r, u], y[0] == r, y'[0] == r*u}, y, {x, 0, 2},
InterpolationOrder -> All]
func2[r_?NumericQ, u_?NumericQ] := (z''[x] == y'[x]*r + u*y[x]*z'[x]) /. sol[r, u]
sol2[r_, u_] := NDSolve[{func2[r, u], z[0] == 1, z'[0] == 0}, z, {x, 0, 2}]

sol2[2, 1]