I have a question for solving t -> Infinity on Mathematica. First, I have a system of ODEs:

feqn1 = {x'[t] == -k1*x[t]*y[t]/h[t] + k2*z[t],
   y'[t] == 3*(-k1*x[t]*y[t]/h[t] + k2*z[t]),
   z'[t] == k1*x[t]*y[t]/h[t] - k2*z[t],
   h'[t] == -3*(-k1*x[t]*y[t]/h[t] + k2*z[t]),
   x[0] == .001, z[0] == 0, y[0] == 0.2, h[0] == .001};

Then I solved the equations by NDSolve.

NDSolve[feqn1, {x, y, z, h}, {t, 0, 1000}];

and for the second part of the problem, I have another set of ODE.

seqn1 = {x2'[t] == -k1*x2[t]*y[t]/h2[t] + k2*z[t],
         h2'[t] == -3*(-k1*x2[t]*y[t]/h2[t] + k2*z[t])}

The initial conditions are h2[0] == 1 at t = Infinity, x2 = 0.001.


  • 2
    $\begingroup$ You can approximate the boundary condition at infinity with a suitably large number, e.g. h2[1/$MachineEpsilon] == 1. Also, see this, among a number of other old MathGroup posts on handling such boundary conditions. $\endgroup$ Nov 8, 2012 at 0:24
  • 4
    $\begingroup$ An alternative is to change variables in your ODE so that infinity becomes 1 say (?) $\endgroup$
    – chris
    Nov 8, 2012 at 8:45
  • $\begingroup$ Thanks J.M & Chris... I tried to look up the document, they used DShoot, but I couldn't find a command on Mathematica for "DShoot". And I also tried set x2[t=1000]= 0.001, but Mathematica only gave me error msg. $\endgroup$
    – DumbleKo
    Nov 8, 2012 at 18:53

1 Answer 1


The comments have made two suggestions - one that we change variables so that $[0,\infty)$ corresponds to a bounded region like $[0,1)$ and the other that we simply solve over a large interval of the form $[0,M]$. Since you've already got InterpolatingFunctions defined over $[0,1000]$ in your system, I suppose that the second suggestion probably makes the most sense. Also, as appealing as a change of variables sounds, in my experience this likely generates a singularlity at the right endpoint making anyway.

Here's your initial system; I grabbed values of $k_1$ and $k_2$ from another post.

k1 = 0.000392;
k2 = 0.00759355499;
feqn1 = {x'[t] == -k1*x[t]*y[t]/h[t] + k2*z[t],
   y'[t] == 3*(-k1*x[t]*y[t]/h[t] + k2*z[t]),
   z'[t] == k1*x[t]*y[t]/h[t] - k2*z[t],
   h'[t] == -3*(-k1*x[t]*y[t]/h[t] + k2*z[t]),
   x[0] == .001, z[0] == 0, y[0] == 0.2, h[0] == .001};
{x,y,z,h} = {x,y,z,h} /. First[

We'll use a shooting method to solve the boundary value problem. That is, we'll write a function shot that solves an initial value problem in terms of a parameter a, with initial conditions at the left endpoint. The idea is to choose the parameter a so that the right hand boundary condition is satisfied. Here's shot:

shot[a_?NumericQ] := Module[{},
  solution = First[NDSolve[
     {x2'[t] == -k1*x2[t]*y[t]/h2[t] + k2*z[t],
      h2'[t] == -3*(-k1*x2[t]*y[t]/h2[t] + k2*z[t]),
      h2[0] == 1, x2[0] == a},
     {h2, x2}, {t, 0, 1000}]];
  {h2Try, x2Try} = {h2, x2} /. solution;
Plot[shot[a], {a, -2, 2}]

enter image description here

The plot indicates that shot is one-to-one and that there's probably a unique solution to your problem. We can find it with FindRoot:

FindRoot[shot[a] == 0.001, {a, 0}]

(* Out: {a -> -0.00479028} *)

Note that the solutions for $x_2$ and $h_2$ are already stored in the variables x2Try and h2Try. We can plot them over various intervals as below. I'm actually plotting $x_2(t)$ and $h_2(t)-1$ here since, otherwise, they look like two flat lines.

 Plot[{x2Try[t], h2Try[t] - 1}, {t, 0, max}],
 {{max, 1}, 0.1, 1000}]

enter image description here

  • $\begingroup$ If need be, one could save on memory in the definition of shot[] by having the start and endpoints in the integration interval coincide: x2[1000] /. First @ NDSolve[(* equations *), {h2, x2}, {t, 1000, 1000}]. $\endgroup$ Nov 8, 2012 at 13:31
  • $\begingroup$ @ Mark McClure: Thanks for the answer, but I'm still trying to figure out "a" because according to your answer, it seems like at t=0, x2=a, which equals to some values. But in my case, at t=0, x2=0 and t=inifinity, x2=0.001 is the boundary conditions. Maybe I misunderstood something here... $\endgroup$
    – DumbleKo
    Nov 8, 2012 at 18:52
  • $\begingroup$ @DumbleKo At no point in your question do I see a specification for x2[0], which is a good thing! You've got a first order system in the two unknown functions h2 and x2 - correct? You've already specified that h2[0]==1. If you additionally specify that x2[0]==0, then your solution is uniquely determined so you can't expect to force the value of x2 at infinity. If you want some particular value for x2 at infinity, then you might be able to get it, if you allow flexibility at zero. That's the approach taken here. $\endgroup$ Nov 8, 2012 at 20:29
  • $\begingroup$ @ Mark McClure: I see. Sorry I was confused by the my question and trying to solve this problems from different approach. But thanks for answering all my questions. They've been very helpful. $\endgroup$
    – DumbleKo
    Nov 8, 2012 at 21:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.