read MarkovWZ:
H:=(z+x)!/(3*x+z+1)!;

t0:=time():

for r from 2 to 6 do
 print(`Using`, H^r, `as a kernel for Zeta(`,r,`) yields the MWZ pair`):
 lu:=Markov(H^r,z,x):
 print(lu):
 print(`The matrix has dimension`, nops(lu[2])):
 print(`The error in the approximation using the 20^th term is`):
 print(evalf(Appxf(H^r,z,x,20)[2]-Zeta(r))):
 print(`The error in the approximation using the 40^th term is`):
 print(evalf(Appxf(H^r,z,x,40)[2]-Zeta(r))):
t1:=time():
print(`This took `, t1-t0, `seconds of CPU time`):
t0:=t1:
od:

quit: