I have the next figure with inflow and outflow rates: figure

And I need to calculate the time for which the pollution will be 50% from the initial value. Assuming that all lakes have the same pollution concentration p initially. The pollution flows from a lake to another in the chain and only inflows not from a lake are clean water.

I wrote the equations for first 3 lakes:

s'[t] + 15/2900 * s[t] == 0

m'[t] + 38/1180 * m[t] == 0

h'[t] + 68/850 * h[t] - 15/2900 * s[t] - 38/1180 * m[t] == 0

And because all the lakes initially have the same pollution concentration p => s(0)==p, m(0)==p, h(0)==p, right?

First equation:

eqS = s'[t] == -15/2900*s[t]
solS = DSolve[{eqS, s[0] == p}, s, t]
s[t_] = s[t] /. First@solS

Second equation:

eqM = m'[t] == -38/1180*m[t]
solM = DSolve[{eqM, m[0] == p}, m, t]
m[t_] = m[t] /. First@solM

Third equation:

eqH = h'[t] == 15/2900*s[t] + 38/1180*m[t] - 68/850*h[t]
solH = DSolve[{eqH, h[0] == p}, h, t]
h[t_] = h[t] /. First@solH

And to find the time for the first and second equation I used:

Solve[s[t] == 0.5*p]
Solve[m[t] == 0.5*p]

And it worked ok, but trying to do the same thing for h[t] takes too much time to compute, and I wonder if my solution is ok.

  • 2
    $\begingroup$ aside to the mathematica content, but I don't think you have the equations right. Why should the huron equation involve the volume of the upstream lakes? $\endgroup$
    – george2079
    Jan 12, 2018 at 16:49
  • $\begingroup$ @george2079 because it says that the pollution goes from a lake to another, the formula is dH/dt = input rate - output rate and for huron the input rate will be the 2 equations for Superior and Michigan. $\endgroup$ Jan 12, 2018 at 22:16
  • $\begingroup$ I get the concept. I just think the third equation should be h'==(15s+38m-68h)/850. That doesn't affect how you solve it of course. $\endgroup$
    – george2079
    Jan 12, 2018 at 22:37

1 Answer 1


The problem is that the function h is too complicated for solving it symbolically with Solve. You can use FindRoot to obtain a numerical solution as follows:

eq = Simplify[h[t]/p == 1/2]
FindRoot[eq, {t, 0}]

{t -> 18.0728}

By the way: You can also solve the system of differential equations at once; this might get handy for more complicated task, e.g. if there were cyclic dependencies between the seas:

ClearAll[s, m, h]
des = {
   s'[t] + 15/2900*s[t] == 0,
   m'[t] + 38/1180*m[t] == 0,
   h'[t] + 68/850*h[t] - 15/2900*s[t] - 38/1180*m[t] == 0
initials = {s[0] == p, m[0] == p, h[0] == p};
{s, m, h} = DSolveValue[Join[des, initials], {s, m, h}, t];

Function[{t}, ( E^(-2 t/25) (7867 + 20615 E^(141 t/2950) + 2115 E^(217 t/2900)) p)/30597]


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.