@bbgodfrey I took your advice above and used Piecewise option. But I did notice that Maybe I had the equation wrong. Please take a look at my added part of my question above.

@bbgodfrey Thank you. I think when s=0, t cannot be less or equal to 2 because of log expression has $log(-2+t)$ in it. Now I want to exclude this $s=0, t\leq 2$ line from my integration (s=0 only). How should I modify the code? I tried to search exclusion option for it, but perhaps I am doing something wrong.

Thank you, but why do I get an error (returns h$1874, with 'Solve::nsmet: This system cannot be solved with the methods available to Solve.') When I plug in a * x * (Log[1 + cx^(1/3) + dx^(2/3)]) as my function to this code? (all I did was multiplying x)

legendreTransform[ Refine[x^2 (Log[1 + x^(1/3) + x^(2/3)]), Assumptions -> {x > 0, Element[x, Reals]}], x, y] doing this will not help either?

@bills Where may I find info. on choosing the branch on Mathematica? is there a built-in function for it?

@bills then is it fine to add change the last part ' Assumptions -> {x>0, Element[x, Reals], Element[x, Reals]}' to solve that problem?

@J.M. I had to examine both cases... sorry

@J.M. This code worked well on the first part. Then I tried to see the behavior when z goes to infinity. Here i my code ord = 3; mzi = M - Exp[-z]*(Sum[a[n] (1/z)^n, {n, 0, ord}] + O[z]^(ord + 1)); pzi = Exp[-z]*(Sum[p[n] (1/z)^n, {n, 0, ord}] + O[z]^(ord + 1)); With[{mzii = mzi, pzii = pzi}, SolveAlways[ D[pzii, z] == -[Epsilon]^2/ z^2 (D[mzii, z]/(4 Piz^2) + pzii) (4*Piz^3*pzii + mzii), z]] Then I want to solve p[n] in terms of a[n]. This gives me error...

Thank you, but what should I do if I want to expand the series at other z (e.g., 1, 2, or infinity), rather than z=0?