Problem of `Integrate` and `DSolve` with inhomogeneous Legendre differential equation

I have the following inhomogeneous Legendre differential equation $(m=n=2)$:

 deq=(6 - 4/(1 - z^2)) h2[z] - 
  2 z Derivative[1][h2][z] + (1 - z^2) Derivative[2][h2][z] == (
  6 (-1 + z + z^2) μ^2)/(M^4 (-1 + z^2)^2) + (

  3 μ^2 Log[(-1 + z)/(1 + z)])/(M^4 (-1 + z))

where the variable $z>1$. I can not get Mathematica to give me the simplest/correct solutions. I think it has to do with those logarithms. I know that the homogeneous solution is

solH=A*LegendreQ[2,2,3,z]+B*LegendreP[2,2,3,z]

I need mathematica's type 3 Legendre functions because for $z>1$ the type 3 functions are on the right branch. With DSolve however I can only get the expanded forms of type 1:

DSolve[deq[[1]] == 0, h2[z], z]

gives back

 {h2[z] -> -3 (-1 + z^2) C[1] + C[2] (((1 - z^2) (5 z - 3 z^3))/(-1 + z^2)^2 + 

      3 (1 - z^2) (-(1/2) Log[1 - z] + 1/2 Log[1 + z]))}

which is not real for $z>1$. Is there a way to tell DSolve to work with $z>1$ or to search for real solutions under that assumption? $Assumptions={z>1} alone does not do the trick.

Apart from that I have a bigger problem with the particular solution: DSolve[deq, h2[z], z][[1]] gives a ridiculous partial solution which is incredibly long and complex. I want a real particular solution; I want a solution that is real for all $z>1$. I know from literature that there is one (plugging that one in proves that it is indeed one):

{h2p[z] -> -3 μ^2/(16 M^4) (6 z^2 + 3 z - 6 - (4 z^2 + 2 z)/(z^2 - 1)) - 3 μ^2/(32 M^4) (3 z^2 - 8 z - 3 - 8/(z^2 - 1)) Log[(z - 1)/(z + 1)] + 3 μ^2/(16 M^4) (z^2 - 1) Log[(z - 1)/(z + 1)]^2}

Again is there a way to tell DSolve to produce a real solution for $z>1$? I also tried constructing that particular solution using the Variation of Parameters method

WLQ223LP223 = 
 Wronskian[{LegendreQ[2, 2, 3, z], LegendreP[2, 2, 3, z]}, z]
-LegendreQ[2, 2, 3, z]*

    Integrate[deq[[2]]*LegendreP[2, 2, 3, z]/WLQ223LP223, z] + 
   LegendreP[2, 2, 3, z]*
    Integrate[deq[[2]]*LegendreQ[2, 2, 3, z]/WLQ223LP223, z] // 
  FunctionExpand // Simplify

Mathematica finds an analytical solutions but they do not reduce to the particular solution I want: the involved integrals become very lengthy and again complex.

Anyone around here who has some experience with such differential equations in Mathematica?

I send the problem to a friend who has Maple and Maple's dsolve finds the "correct"/simple partial solution, when one prepends assume(z > 1). It does so using the Variation of Parameters method and the solutions of the homogeneous equation.

Answer

Personally, I prefer @corey979's approach (+1), but with the following tweaks.

solExpr = 
  Assuming[z > 1 && M ∈ Reals && μ ∈ Reals,
   DSolveValue[deq, h2[z], z] /. {Log[-1 + z^2] -> Log[-1 + z] + Log[1 + z], 
        Log[(-1 + z)/(1 + z)] -> Log[-1 + z] - Log[1 + z]} // Re // 
     ComplexExpand[#, TargetFunctions -> {Re, Im}] & // Simplify
   ];

sol = h2 -> Function @@ {z, solExpr}

deq /. sol // Simplify[#, z > 1 && M ∈ Reals && μ ∈ Reals] &
(*  True  *)

Plot[h2[z] /. sol /. {M -> 1, μ -> 1, C[1] -> -1, C[2] -> 1},
 {z, 1, 5}]

The main tweaks are including all assumptions to make all parameters real in both DSolve and Simplify and the use of TargetFunctions. It's also easier to simplify an algebraic expression; hence DSolveValue[]. The rules

{Log[-1 + z^2] -> Log[-1 + z] + Log[1 + z], 
 Log[(-1 + z)/(1 + z)] -> Log[-1 + z] - Log[1 + z]}

help with simplification; they are not important and may be omitted.

The main reason I prefer this is that the real and imaginary parts are independently solutions of a real ODE, a key point made by corey979. It just seems a simpler approach, when it simplifies down nicely. In this case it does, and we get a general solution that spans the solution space.