# NDSolve: obtain time difference when different solutions get same distance to equilibrium value

I have a numerically solved a differential equation which has an equilibrium value and a starting point, both included as a parameter. I want to plot the difference between the times it takes to come into a certain proximity of the equilibrium value starting from +/- starting parameter. I'll show you the code and hopefully things get more clear:

bsp[c_, d_] := NDSolve[{x'[t] == -x[t] + d, x[0] == c}, x, {t, 0, 30}]
(* defines the differential eqn *)
bspdiff[d_] := Module[{d1 = d, func, func2, diff, par, x, t, t1},
func := x /. First[bsp[par, d1]];
func2 := x /. First[bsp[-par, d1]];
Plot[
Abs[(t /. FindRoot[func[t] == 1.01*d1, {t, 2}]) - (t1 /. FindRoot[func2[t1] == 0.99*d1,
{t1, 2}])],
{par, 1, 5}]]


As you see, I define two functions, one for the positive, one for the negative value of the starting point. Then I use FindRoot to obtain the times at which they are within 1% of the equilibrium value and deduct those two. I want to plot this result as a function of the starting point. What I get upon executing this module is an empty plot and:

FindRoot::nlnum: The function value {-2.02+x\$70015[2.]} is not a list of numbers with dimensions {1} at {t\$70015} = {2.}. >>

ReplaceAll::reps: {FindRoot[func\$70015[t\$70015]==1.01 2,{t\$70015,2}]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. >> FindRoot::nlnum: The function value {-1.98+x\$70015[2.]} is not a list of numbers with dimensions {1} at {t1\$70015} = {2.}. >> ReplaceAll::reps: {FindRoot[func2\$70015[t1\$70015]==0.99 2,{t1\$70015,2}]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. >>

FindRoot::nlnum: The function value {-1.98+x\$70015[2.]} is not a list of numbers with dimensions {1} at {t1\$70015} = {2.}. >>

General::stop: Further output of FindRoot::nlnum will be suppressed during this calculation. >>

ReplaceAll::reps: {FindRoot[func2\$70015[t1$70015]==0.99 2.,{t1\\$70015,2.}]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. >>

General::stop: Further output of ReplaceAll::reps will be suppressed during this calculation. >>

When I however execute this manually for single values, it works:

func := x /. First[bsp[2, 1]];
func2 := x /. First[bsp[-2, 1]];
Abs[(t /. FindRoot[func[t] == 1.01, {t, 2}]) - (t1 /.
FindRoot[func2[t1] == 0.99, {t1, 2}])]


yields a numerical output. What I do not understand - where is my mistake/why does it not work? Or is there a better way to do this? I would be very happy about any help and suggestions.

Here's how I would do this, using WhenEvent to find the times you're looking for while solving the differential equation and Sow and Reap to extract them.

Clear[bsp]
bsp[c_ /; c > 1, d_] := With[
{t1 = First @@ Last@Reap@NDSolve[{
x'[t] == -x[t] + d
, x[0] == c
, WhenEvent[x[t] - 1.01 d == 0, {Sow[t], "StopIntegration"}]
}
, x
, {t, 0, 5000}]
, t2 = First @@ Last@Reap@NDSolve[{
x'[t] == -x[t] + d
, x[0] == -c
, WhenEvent[x[t] - 0.99 d == 0, {Sow[t], "StopIntegration"}]
}
, x
, {t, 0, 5000}]}
, Abs[t1 - t2]
]


Then, just do

Plot[bsp[c, 1], {c, 1.01, 5}]


resulting in

## Pre-Version 9

WhenEvent was introduced in MMA Ver. 9. In previous version, we can use the NDSolve Method called "EventLocator". Here is the adapted code:

Needs["DifferentialEquationsNDSolveUtilities"];

Clear[bsp]
bsp[c_ /; c > 1, d_] := With[{
t1 = First @@ Last@Reap@NDSolve[{x'[t] == -x[t] + d, x[0] == c}
, x
, {t, 0, 10000}
, Method -> {"EventLocator"
, "Event" -> (x[t] == 1.01 d)
, "EventAction" :> (Sow[t]; "StopIntegration")}
]
, t2 = First @@ Last@Reap@NDSolve[{x'[t] == -x[t] + d, x[0] == -c}
, x
, {t, 0, 10000}
, Method -> {"EventLocator"
, "Event" -> (x[t] == 0.99 d)
, "EventAction" :> (Sow[t]; "StopIntegration")}
]
, Abs[t1 - t2]
]


## Analytic version

Of course, in this particular case, there is an analytic solution to the problem.

sols = x[t] /. First@DSolve[{x'[t] == -x[t] + d, x[0] == c}, x[t], t]
times = t /. {Solve[sols == 1.01 d, t], Solve[sols == 0.99 d /. c -> -c, t]} // Flatten // Quiet
Plot[
Evaluate[{sols, sols /. c -> -c, 1.01 d, 0.99 d} /. {d -> 1, c -> 2}]
, {t, 1, 10}
, GridLines -> {times /. {d -> 1, c -> 2}, None}]


results in

Then, we can define

Clear[bsp]
bsp[c_ /; c > 1, d_] = Abs@*Subtract @@ times;


or, pre MMA-10

bsp[c_ /; c > 1, d_] = Abs~Composition~Subtract @@ times;

• I am very sorry - I did not include that I am working with an old version of Mathematica (8) as it did not seem relevant to me. I am however very thankful for your answer and might just need to obtain a newer version then ... Sep 1, 2015 at 15:11
• @Inari. A similar method exists for older versions! I will write that up, although I have to remember how it works... Sep 1, 2015 at 15:15
• @Inari. I believe the second answer will work for you. See updated post. You call bsp[c, d] in the same way as before. Sep 1, 2015 at 15:36
• I cannot thank you enough. It's wonderful and it works perfectly well. Sep 1, 2015 at 15:51
• @Inari. Thank you for the accept! But it is a good idea to wait for at least a day before accepting just in case someone else wants to post an alternative (and possibly better!) answer. When people see a question with an accepted answer, they less likely to post additional answers; questions with unaccepted answers get more answers! That said, I believe that this is the best way to go about solving this problem. Sep 1, 2015 at 16:07