If not for the positive integer restriction, (x,y,z) = (0,t,t) or (t,0,t) would work, but x,y,z must be non zero, so scratch that one.

After some algebra, one parametric solution is:

(x,y,z) = (t-1, t, t^2)

(x,y,z) = (t, t-1, t^2)  (equivalent since x and y are interchangeable)

Proof:

Consider the substitution of y = (x+1)

(x^2+x+1)*((x+1)^2+(x+1)+1) = 

(x^2+x+1)*(x^2 + 2x + 1+x+2) = 

(x^2 + x + 1)*(x^2 + 3x + 3) = 

x^4 + 4x^3 + 7x^2 + 6x + 3 = 

(x^4 + 4x^3 + 6x^2 + 4x + 1) + (x^2 + 2x + 1) + 1

(x+1)^4 + (x+1)^2 + 1

letting z = (x+1)^2 yields z^2+z+1

But wait, there's more.

Using a program, (setting WLOG x <= y), I found several other solutions but I have not been able to put them into a parametric form.  Note that I did not find any other solutions other than the above parametric formula for any x larger than 11, but that might not hold if larger y and z values were allowed.

x y z

1 2 4

1 9 16

1 35 61

1 132 229

1 494 856

2 3 9

2 11 30

2 25 67

2 56 149

2 183 485

2 406 1075

3 4 16

3 49 178

3 87 315

3 305 1101

4 5 25

4 26 121

4 81 373

4 390 1789

5 6 36

5 84 470

5 121 676

6 7 49

6 169 1111

6 230 1511

7 8 64

7 74 562

7 225 1702

8 9 81

9 10 100

9 55 529

9 68 653

10 11 121

10 146 1543

11 12 144

11 129 1493

12 13 169

13 14 196

14 15 225

 t-1  t   t^2  afterwards

Putting in the y values when x=1:   2,9,35,132,494  gives oeis A001571

a(n) = 4*a(n-1) - a(n-2) + 1, with a(0) = 0, a(1) = 2

While putting in the z values when x=1:  4,16,61,229,856 gives oeis A133161

Indices of the triangular numbers which are also centered triangular number.

This is not a complete solution.  For one thing, allowing larger values for y and z will find more solutions.  For another, I don't have a parametric version for these other solutions.