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}]
