How to recover y-coordinates when using XZ Montgomery curve

Question

I am using Montgomery ladder with Montgomery curve $by^2=x^3+ax^2+x$ using XZ coordinates and I recovered the $X$ value using $X3=X1/Z1$, but I don't know how to recover the $Y$ coordinates.

for Double and add ladder I am using this:

      A = X2+Z2
      AA = A2
      B = X2-Z2
      BB = B2
      E = AA-BB
      C = X3+Z3
      D = X3-Z3
      DA = D*A
      CB = C*B
      X5 = Z1*(DA+CB)2
      Z5 = X1*(DA-CB)2
      X4 = AA*BB
      Z4 = E*(BB+a24*E)

I tried this way :

x3=2;
y3 = mod(mod((x3.^3 + mod(a*x3.^2,p)+x3),p) * mod(modinvr(b,p),p),p);

for y = 0:22
    x = mod(y^2, 23);
    if x == y3
        fprintf("y = %d\n", y);// here I got two values of y 8 and 15
    end
end

here I got two values of y 8 and 15 both are correct points on the curve but in my case I want to choose 8 because the affine scalar point is (2,8) I have another point on the curve (2,15) but not in my scalar point! so that's why I need to select 8 instead of 15.

values used: a=6 b=7 base point (8,1) and these the points on the curve:

     0     0
     2     8
     2    15
     3     9
     3    14
     7     0
     8     1
     8    22
     9     2
     9    21
    10     0
    11     7
    11    16
    12     2
    12    21
    17     3
    17    20
    18     4
    18    19
    19     2
    19    21
    21     5
    21    18

the resource is from this link: Montgomery curve XZ coordinates