restart;
convert(harmonic(b), Sum):

    Sum(b/_k1/(_k1+b),_k1 = 1 .. infinity)

value(%);

               gamma + Psi(1 + b)

T1 := expand(value(convert(harmonic(z),Sum)) = harmonic(z));

                    1
     T1 := Psi(z) + - + gamma = harmonic(z)
                    z

Q1 := Sum(sum(1/k, k=1..b), b=1..m);

       Sum(gamma+Psi(1+b),b = 1 .. m)

simplify(subs(b=z,  combine(expand(Q1))), {T1});

       Sum(harmonic(z),z = 1 .. m)

value(%);

      harmonic(m + 1) (m + 1) - m - 1

T2 := isolate(T1, Psi(z));

                                          1
     T2 := Psi(z) = harmonic(z) - gamma - -
                                          z

temp := sum(sum(1/k, k=1..b), b=1..m):

Q2 := simplify(expand(temp),size);

         / 2                               2
         m  + m/ Psi(m) + 1 + (gamma - 1) m  + (gamma + 1) m
   Q2 := ----------------------------------------------------
                                  m

R := subs(Psi(z)=Psi(z::anything), T2):

normal(applyrule(R, Q2));

        harmonic(m) m + harmonic(m) - m

normal(subs(subs(z=m, T2), Q2));

        harmonic(m) m + harmonic(m) - m