from sympy import *
init_printing()
p, r, c, p, y, lam, f = symbols('p r c p y lambda f')
constant1 = exp(-lam) / (1 - exp(-lam))
constant1
constant2 = constant1 * (exp(lam) - 1 - lam)
constant2
eta = 1 - constant2 + constant1 * (exp(lam * (1-f)) - 1 - lam * (1 - f))
eta
etaOfR = limit(eta, f, 1)
etaOfR
expression1 = Eq(y / p * (p - c) * eta,
y / r * (r - c) * etaOfR)
expression1
sol = solve(expression1, f)
sol
sol[0].simplify()
This is the result from Mathematica
$f\to \frac{r (c-p) W\left(-e^{\frac{r (c (\lambda -1)+p)-c \lambda p}{r (c-p)}}\right)+c (\lambda p+r)-(\lambda +1) p r}{\lambda r (c-p)}$