title: Can it be that `sympy` is much, much slower than Mathematica?
body: <p>I'm reproducing Mathematica results using <code>Sympy</code>, and I'm new to the latter, so I might be doing things wrong. However, I noticed that some stuff that took a minute at max using Mathematica is just taking forever (read: did not finish after I started it an hour ago) in sympy. That applies both to <code>Simplify()</code>, and <code>solve()</code>. Am I doing something wrong, or is that really the case?</p>
  
  <p>I'll attach my <code>solve()</code> case:</p>
  
  <pre><code>import sympy as sp
  from sympy import init_printing
  init_printing()
  
  p, r, c, p, y, Lambda = sp.symbols('p r c p y Lambda')
  
  F = sp.Symbol('F')
  eta1 = lambda p: 1/(1-sp.exp(-Lambda) * sp.exp(-Lambda)*(sp.exp(Lambda) - 1 - Lambda))
  eta2 = lambda p: 1/(1-sp.exp(-Lambda)) * sp.exp(-Lambda)/(1-F) * (sp.exp(Lambda*(1- F)) - 1 - Lambda*(1-F))
  
  eta = lambda p: 1 - eta1(p) + eta2(p)
  etaOfR = sp.limit(eta(p), F,  1)
  
  S = lambda p: eta(p)*y/p*(p-c)
  SOfR = etaOfR*y/r*(r-c)
  sp.solve(S(p)-SOfR, F)
  </code></pre>
  
  <p><em>The corresponding Mathematica code</em>:</p>
  
  <pre><code>ClearAll[r, p, lambda, a, A, c, eta, f, y, constant1, constant2, eta, \
  etaOfR]
  constant1[lambda_] := Exp[-lambda]/(1 - Exp[-lambda]);
  constant2[lambda_] := constant1[lambda]*(Exp[lambda] - 1 - lambda);
  eta[lambda_, f_] := 
    1 - constant2[lambda] + 
     constant1[lambda]*(Exp[lambda*(1 - f)] - 1 - lambda*(1 - f)) ;
  etaOfR[lambda_] := Limit[eta[lambda, f], f -&gt; 1];
  expression1[lambda_, f_] := 
    y/p (p - c) eta[lambda, f] == y/r (r - c) etaOfR[lambda];
  
  Solve[expression1[lambda, f], f] // FullSimplify
  </code></pre>
  
  <p>Output:</p>
  
  <pre><code>{{f -&gt; (-(1 + lambda) p r + 
      c (lambda p + r) + (c - 
         p) r ProductLog[-E^(((-c lambda p + (c (-1 + lambda) + 
             p) r)/((c - p) r)))])/(lambda (c - p) r)}}
  </code></pre>