# How to solve in principle such coupled ODEs

One example for coupled ODEs from the Mathemtica help center reads

a = NDSolve[{x'[t] == -3 (x[t] - y[t]),
y'[t] == -x[t] z[t] + 26.5 x[t] - y[t],
z'[t] == x[t] y[t] - z[t], x == z == 0,
y == 1}, {x, y, z}, {t, 0, 5}]


Now we can for example extend this expression by replacing

x'[t] == -3 (x[t] - y[t]) -> x['t] == -3* (x[t] - y[t])*int[t]


with

f[t_] := Sin[t^2]
int[t_?NumericQ] := NIntegrate[Exp[-f[t]]*v^2, {v, 0, 1}]


This works indeed fine, but now I would like to extend the integral expression further by making the integrand dependent on the current solution of z[t] i.e.

int[t_?NumericQ] := NIntegrate[Exp[-f[t]*Evaluate[z[t]/.a]*v^2], {v, 0, 1}]


But this expression returns a bunch of errors. Does anyone know how to make this work?

• But neither z nor f are dependent on v, no? Jul 25 '16 at 13:41
• Right, they are both independent
– stef
Jul 25 '16 at 13:59
• Then you can factor them out of the integral, no? Jul 25 '16 at 14:21
• Thats not trivial if one has e.g. Sinh[f[t]*z[t]*v] as integrand
– stef
Jul 25 '16 at 14:55
• You probably should have used that example, then. :) Jul 25 '16 at 14:58

I think substituting a in the definition of int, which is called by a is trouble, and what you really need is to make sure int isn't evaluated when z[t] isn't a number. The following seems to work, albeit not very quickly:

Clear[a];

f[t_] := Sin[t^2];
int[t_?NumericQ, z_?NumericQ] := NIntegrate[Exp[-f[t]*z*v^2], {v, 0, 1}];

a = NDSolve[{
x'[t] == -3*(x[t] - y[t])*int[t, z[t]],
y'[t] == -x[t] z[t] + 26.5 x[t] - y[t],
z'[t] == x[t] y[t] - z[t],
x == z == 0, y == 1},
{x, y, z}, {t, 0, 5}]

• That's what I was looking for, works fine! However, adding to the above definition of int[t,z] the new def. int2[t_?NumericQ, z_?NumericQ] := h /. NDSolve[{h'[v] == Exp[-f[t]*z*v^2], h == 0}, h, {v, 0, 1}] and calling int2[t,z[t]] in NDSolve produces fatal errors, even though Plot[{int[t, 1], int2[t, 1]}, {t, 0, 5}] shows the same result. Do you know how to do the replacement right?
– stef
Jul 27 '16 at 9:14
• int2[2, 1] gives {1.32159} not 1.32159 -- you have to strip off a level with []. That is, int2[t_?NumericQ, z_?NumericQ] := h /. NDSolve[{h'[v] == Exp[-f[t]*z*v^2], h == 0}, h, {v, 0, 1}][] works Jul 27 '16 at 18:38