I found it interesting that with some algebraic manipulations the system of the two differential equations can be decoupled to give an expression for y as a function of f.
yy = (a (-f + f0))/b + y0 + Log[f] - Log[f0]
(Note : you should run the code given above to have all symbols available).
With appropriate modifications this is indeed the same as fy given above, or better the difference is reasonably well near to zero:
Plot[fy[0] - .2/.05 (ff[t] - ff[0]) + Log[ff[t]/ff[0]] - fy[t], {t, 0, 240},
PlotRange -> All]
Now you can insert this y in the 2nd equation and obtain an equation for f
fp1 = - b yy f // Simplify
dgl = fF'[t] == (fp1 /. {f -> fF[t], f0 -> ff[0], y0 -> fy[0], a -> .2, b -> .05})
Unfortunately this cannot be integrated, but NDSolve does the job
Clear[fF]
fF = fF /. Flatten[NDSolve[ {dgl , fF[0] == .9999}, fF, {t, 0, 200}]]
fF is equivalent to ff
Plot[{ff[tt], fF[tt]}, {tt, 0, 150},
PlotStyle -> {Blue, {Thick, Dashing[{.05, .05}], Red}}]
and we can insert it into the expression for y to obtain
fY = (yy /. {a -> .2, b -> .05, y0 -> fy[0], f0 -> ff[0], f -> fF[t]})
which as well is equivalent to fy obtained earlier :
Plot[{fY, fy[t]}, {t, 0, 180},
PlotStyle -> {Blue, {Thick, Dashing[{.05, .05}], Red}}]
If there is the need to operate with some more explicit functions one could fit the derivative of ff (or fF) to some appropriate functions. Here I chose a Gaussian and a Lorentzian
data = Table[{t, -D[ff[x], x] /. x -> t}, {t, 50, 150}];
lsg = FindFit[
data,
p1 Exp[-p2 (t - p3)^2] + p4/(1 + p5 (t - p6)^2),
{{p1, .015}, {p2, .005}, {p3, 84}, {p4, .015}, {p5, .005}, {p6,
84.5}},
t]
giving
empfp = p1 Exp[-p2 (t - p3)^2] + p4/(1 + p5 (t - p6)^2) /. lsg
which looks quite reasonable :
Plot[empfp, {t, 40, 120}, Epilog -> Point /@ data]
This can be integrated to give an approximation for ff :
SetOptions[Integrate,GenerateConditions->False];
empff = Integrate[-empfp, {t, 0, x}] + 0.9999
Plot[empff, {x, 0, 140},
Epilog -> Point /@ Table[{t, ff[t]}, {t, 20, 140, 5}]]
Insert this in the first differential equation and integrate to get fy
empfy = Exp[ Integrate[(a empff - b /. a -> .2 /. b -> .05), {x, 0, t}]]
The leading factor must be adjusted a bit to give a nice result
Plot[1.282 10^-5 empfy, {t, 0, 250},
Epilog -> Point /@ Table[{t, fy[t]}, {t, 0, 200, 2}],
PlotRange -> All]
Attachments: