No algebraic solution exists, at least not that it could find. I'm not at all surprised at that, Th is inside and outside of the acosh, and the presence of sqrts in there will cause more trouble. The solution you got is an indication that no solution was found, just expressed in a confusing way.
First, skip the part where you set those equations to some value anyway. You need to be able to plot things. And you cannot plot an equation expressed as f(x)==0. You just want to plot f(x), and see if it ever crosses zero.
We can simplify the problem. ls really is not needed as a parameter. Subtract equations 1 and 2, and 2 and 3. We want to find a solution for the 2 variable system of equations. Does one exist?
vpa(isolate(equ2-equ3,dm),5)
ans =
dm == 2.0844e-13*Th - 6.5234e-45*(1.021e+63*Th^2 + 2.4338e+72*Th + 3.9302e+71)^(1/2) + 0.00024844
So I solved for dm, as a function of Th.
Do you see anything interesting? I do. Those immense powers of 10 tell me you have a truly serious numerical problem.
For which values of Th does a real solution exist?
format long g
roots([1.021e+63 2.4338e+72 3.9302e+71])
So, in the interval [-2383741429.80913, -0.161484098950869], the inside of the sqrt is negative. dm will be complex in that case. Therefore, if dm must be real and positive, we need for Th to be at least -0.16. You also said that Th must be always positive.
fplot(@(Th) 2.0844e-13*Th - 6.5234e-45*sqrt(1.021e+63*Th^2 + 2.4338e+72*Th + 3.9302e+71) + 0.00024844,[0,1e10])
It looks like dm goes negative when Th exceeds roughly 2e11. So that is the uper limit on Th.
I was never able to find ANY solutions for the difference of equations 1 and 2.
You either have a nasty numerical problem, where the variables vary by many orders of magnitude, or possibly you wrote down the wrong expressions.
