# Solving a system of equations is taking too long

I am trying to find a solution to a fairly simple maximization problem under various modifications, and Mathematica keeps getting stuck on the last one (which seems rather similar to the first, so I'm not sure what the problem is)... Here's the code to my problem:

pp = α*p + (1 - α)

ulpm = pp*Log[clpm] + Log[dlpm] + Log[llpm]
bclpm = wl*(1 - llpm) (1 - tpm) + Tpm - clpm - (1 + τpm)*dlpm
eqV = V == ulpm + λpm*bclpm
foc37 = D[eqV, clpm]
foc38 = D[eqV, dlpm]
foc39 = D[eqV, llpm]
foc40 = D[eqV, λpm]
sols = Solve[{foc37, foc38, foc39, foc40}, {clpm, dlpm,
llpm, λpm}]; {clpm, dlpm, llpm, λpm} /. sols[[1]]
Clear[clpm, dlpm, llpm, λpm]
Set @@@ sols[[1]]

uhpm = pp*Log[chpm] + Log[dhpm] + Log[lhpm]
bchpm = wh*(1 - lhpm) (1 - tpm) + Tpm - chpm - (1 + τpm)*dhpm
eqW = W == uhpm + γpm*bchpm
foc41 = D[eqW, chpm]
foc42 = D[eqW, dhpm]
foc43 = D[eqW, lhpm]
foc44 = D[eqW, γpm]
sols = Solve[{foc41, foc42, foc43, foc44}, {chpm, dhpm,
lhpm, γpm}]; {chpm, dhpm, lhpm, γpm} /. sols[[1]]
Clear[chpm, dhpm, lhpm, γpm]
Set @@@ sols[[1]]

Wpm = pp*Log[clpm] + pp*Log[chpm] + 2*(1 - p)*Log[(clpm + chpm)/2] +
Log[dlpm] + Log[dhpm] + Log[llpm] + Log[lhpm]
rcpm = tpm*(wh*(1 - lhpm) + wl*(1 - llpm)) + τpm*(dhpm + dlpm) -
2*Tpm
eqX = X == Wpm + μpm*rcpm
foc45 = D[eqX, tpm]
foc46 = D[eqX, τpm]
foc47 = D[eqX, Tpm]
foc48 = D[eqX, μpm]
sols = Solve[{foc45, foc46, foc47, foc48}, {tpm, τpm,
Tpm, μpm}]; {tpm, τpm, Tpm, μpm} /. sols[[1]]
Clear[tpm, τpm, Tpm, μpm]
Set @@@ sols[[1]]


The first blocks work fine but the last one has been running for over an hour and has not yet yielded any solution? Any ideas how I can help simplify the equations to get a solution? By the way, a somewhat simpler version of the problem runs smoothly. Here's the code to that problem:

ulpn = p*Log[clpn] + Log[dlpn] + Log[llpn]
bclpn = wl*(1 - llpn) (1 - tpn) + Tpn - clpn - (1 + τpn)*dlpn
eqO = O == ulpn + λpn*bclpn
foc13 = D[eqO, clpn]
foc14 = D[eqO, dlpn]
foc15 = D[eqO, llpn]
foc16 = D[eqO, λpn]
sols = Solve[{foc13, foc14, foc15, foc16}, {clpn, dlpn,
llpn, λpn}]; {clpn, dlpn, llpn, λpn} /. sols[[1]]
Clear[clpn, dlpn, llpn, λpn]
Set @@@ sols[[1]]

uhpn = p*Log[chpn] + Log[dhpn] + Log[lhpn]
bchpn = wh*(1 - lhpn) (1 - tpn) + Tpn - chpn - (1 + τpn)*dhpn
eqP = P == uhpn + γpn*bchpn
foc17 = D[eqP, chpn]
foc18 = D[eqP, dhpn]
foc19 = D[eqP, lhpn]
foc20 = D[eqP, γpn]
sols = Solve[{foc17, foc18, foc19, foc20}, {chpn, dhpn,
lhpn, γpn}]; {chpn, dhpn, lhpn, γpn} /. sols[[1]]
Clear[chpn, dhpn, lhpn, γpn]
Set @@@ sols[[1]]

Wpn = p*Log[clpn] + p*Log[chpn] + 2*(1 - p)*Log[(clpn + chpn)/2] +
Log[dlpn] + Log[dhpn] + Log[llpn] + Log[lhpn]
rcpn = tpn*(wh*(1 - lhpn) + wl*(1 - llpn)) + τpn*(dhpn + dlpn) -
2*Tpn
eqQ = Q == Wpn + μpn*rcpn
foc21 = D[eqQ, tpn]
foc22 = D[eqQ, τpn]
foc23 = D[eqQ, Tpn]
foc24 = D[eqQ, μpn]
sols = Solve[{foc21, foc22, foc23, foc24}, {tpn, τpn,
Tpn, μpn}]; {tpn, τpn, Tpn, μpn} /. sols[[1]]
Clear[tpn, τpn, Tpn, μpn]
Set @@@ sols[[1]]


As a general tip, I would advice to specify assumptions for the symbols you're using. Mathematica tends to give solutions over complex numbers, so if you use logarithms and square roots and such, Mathematica tends to spend a lot of time figuring out what to do with negative/complex values of your variable. Unless you're interested in those complex values, it's better to tell Mathematica to only look in the regions your interest in by using Assuming and \$Assumptions.
Another thing to keep in mind is that it's good to do intermediate simplification using Simplify, FullSimplify etc. This can cut down the complexity of intermediate results a lot, since Solve doesn't simplify by itself.
Finally, Reduce is a better general-purpose tool for equation solving. From the documentation:
Reduce is the thing you want to use whenever Solve seems to get stuck. Reduce has the added bonus that it tells you what assumptions accompany each solution, so you can appreciate why Mathematica spends much time on your solution (you may not be giving it some trivial assumption it needs to simplify the result, for example). And for maximization problems, you may also want to take a look at Maximize since Maximize also checks that your solution is actually a maximum (and not, e.g. just a saddle point) and can deal with singular points as well.
• To include assumptions in your calculation, simply wrap Solve or Reduce in Assuming (got that wrong in my initial post, I corrected that). For example, compare: Assuming[y > 0 && a > 0, Simplify[Solve[a Log[x] == y, x]]] against Solve[a Log[x] == y, x]. Wrap Simplify or FullSimplify (careful: FullSimplify can take a while to evaluate) around every intermediate result you get to make sure that all of the clutter gets cut down. – Sjoerd Smit May 2 '17 at 9:42