x²z + y²z + 4xy = 40,

x² + y² + xyz = 20

Mulitply the first equation by z and the second by 4 and subtract:

x²(z²-4) + y²(z²-4) = 40z-80 
x²(z+2)(z-2) + y²(z+2)(z-2) = 40(z-2)
x²(z+2) + y²(z+2) = 40
(x²+y²)(z+2) = 40
As the first factor must be positivewe only need to check the positive factors of 40.
(x²+y²)=40 (z+2)=1 yields (2,6,-1) which doesn't work in the original system yielding 8 and 28 respectively,
(x²+y²)=20 (z+2)=2 yields (2,4,0) which doesn't work in the original system yielding 32 and 20 respectively,
(x²+y²)=10 (z+2)=4 yields (1,3,2) which doesn't work in the original system yielding 32 and 16 respectively,
(x²+y²)=5 (z+2)=8 yields (2,1,6) which doesn't work in the original system yielding 38 and 17 respectively,
(x²+y²)=4 (z+2)=10 yields (2,0,8) which doesn't work in the original system yielding 32 and 4 respectively,
(x²+y²)=1 (z+2)=40 yields (0,1,38) which doesn't work in the original system yielding 31 and 1 respectively,

but (x²+y²)=2 (z+2)=18 yields <b>(1,1,18)</b> which works.

I think this is the only solution, but looking at the original problem it still seems that a negative value for z could also give solutions.