# System of ODEs - NDSolve issues

I am a self-taught beginner trying to use Mathematica for the first time. If you wouldn't mind, I would like to ask for help with the code I am working on as I keep running into multiple issues when using NDSolve to solve for a problem where given four ODEs and some terminal conditions, I am asked for the time (tf) at which the terminal conditions should be reached.

The equations are somewhat hypothetical for tractability of our real equations that deal with water, so there's not much of a meaningful story behind them.

The issues I face are: Depending on the choice on the FindRoot option for the starting point (t0 + x0), I end up with NDSolve::ndsz, NDSolve::ndsv:, NDSolve::mxst, etc.

Also, I believe that the equations are stiff, and should have a singularity at W=0.

The main problem, however, is that despite all attempts, I still couldn't get any value for tf to solve the problem

With these, I was wondering if anyone could help. Your input would be greatly appreciated!

(*Parameters and Functions*)
t0 = 1991; r = .03; pb = 10; W0 = 5000; s = .05; w = 5000; cq = 3; alpha = 1; g = .02; eta = .3; cd = .947

l = 0.00005*W[t]^2
lprimew = D[l, W[t]]
evap = -.00004 W[t]^2
evapprimew = D[evap, W[t]]
ct = 10*S[t]
ctprime = D[ct, S[t]]
damage = 10*S[t]^2
damageprime = D[damage, S[t]]
q = alpha*E^(g (t - t0)) (p[t] + cd)^(-eta)

(*Terminal W, S, Mu*)
Wprox = W[t] /. Solve[w - l - evap == 0, W[t]]
Wmax = First[Select[W[t] /. Solve[w - l - evap == 0, W[t]], # > 0 &]]

Optans = Select[{\[Mu][t], W[t], S[t]} /. Solve[{pb - cq - ct == -1/(evapprimew + lprimew + r)*((S[t]*damageprime/W[t]) + ((s*w - S[t]*l)/W[t])*ctprime + (\[Mu][t]/(W[t])^2)*(2*s*w - S[t]*(w - evap + evapprimew*W[t]))), \[Mu][t] == -(ctprime*(w - evap - l) + damageprime)/(r + (w - evap)/W[t]), S[t] == s*w/(w - evap - 2*l)}, {\[Mu][t], W[t], S[t]}], Re[#[[1]]] < 0 && Wmax >=  Re[#[[2]]] > 0 && 1 >= Re[#[[3]]] > 0 &]
Wopt = Optans[[1, 2]]
Sopt = Optans[[1, 3]]
\[Mu]opt = Optans[[1, 1]]

(*Differential Equations*)
peqn = p'[t] == (p[t] - cq - ct)*(evapprimew + lprimew + r) + (((S[t]*damageprime/W[t]) + ((s*w - S[t]*l)/W[t])*ctprime + (\[Mu][t]/(W[t])^2)*(2*s*w - S[t]*(w - evap + evapprimew*W[t]))));
Weqn = W'[t] == w - q - evap - l;
Seqn = S'[t] == (s*w - S[t]*(q + l))/W[t];
mueqn = \[Mu]'[t] == r*\[Mu][t] + ctprime*q + damageprime + \[Mu][t]*(q + l)/W[t];

(*Numerical Solution*)
system[x_] := {peqn, Weqn, Seqn, mueqn, p[x] == pb, W[x] == Wopt, S[x] == Sopt, \[Mu][x] == \[Mu]opt}
solution[tf_] := NDSolve[system[tf], {p[t], W[t], S[t], \[Mu][t]}, {t, t0, tf}, Method -> {"StiffnessSwitching", "NonstiffTest" -> False}, AccuracyGoal -> 4, PrecisionGoal -> 4]
Winitial[tf_?NumberQ] := First[(W[t] /. solution[tf]) /. t -> t0]
tfopt = tf /. FindRoot[Winitial[tf] == W0, {tf, t0 + 8}]

Result = solution[tfopt];
Plot[Evaluate[W[t] /. Result], {t, t0, tfopt}]

• Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) 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, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. Dec 30, 2015 at 18:20

A Plot of Winitial

Plot[Winitial[tf], {tf, t0, t0 + 2}]


shows that it is real over only a narrow range of tf, and FindRoot in combination with NDSolve has difficulty outside that range. So, try

tfopt = tf /. FindRoot[Winitial[tf] == W0, {tf, 1991.5}]
(* 1991.73 *)


which is, I believe, the correct answer.

• Thank you so much!! I tried initially higher values as starting points (such as t0 + 9), based on my computations on optimal q, and then moved in multiples of 5 (tried 4). I didn't expect it to be such a small number! I just want to ask how you thought of such a small starting number to try for the FindRoot command. Anyway, thank you again!!! Dec 31, 2015 at 1:54
• @Anonymous I read a guess off the plot. Dec 31, 2015 at 4:00