I will offer a "two stage" manner to tackle this. The first one requires some knowledge of four vectors and familiarity with properties of Minkowski space in the context of SR. The second one I will only outline, one can make it much more formal, I will only describe its intuitive aspect to complete the first more "formal" argument.
If we let the four velocities of the particles, in the lab frame be:
$$ (v^{\mu}_1)^{\text{lab}} = \gamma^{\text{lab}}_{v_1}(1, \vec{v}^{\text{ lab}}_1) $$
$$ (v^{\mu}_2)^{\text{lab}} = \gamma^{\text{lab}}_{v_2}(1, \vec{v}^{\text{ lab}}_2) $$
then we can use the fact that the Minkowskian dot product $$\eta_{\mu\nu}v^{\mu}_1 v^{\nu}_2\tag{1}$$ is invariant across all Lorentz frames. Therefore in the frame where particle $1$ is at rest we have:
$$ (v^{\mu}_1)^{\text{[1]}} = (1, \vec{0}) $$
$$ (v^{\mu}_2)^{\text{[1]}} = \gamma_{v_2}^{[1]}(1, \vec{v}_2^{[1]}) $$
and similarly in the Lorentz frame where particle $2$ is at rest we have:
$$ (v^{\mu}_1)^{\text{[2]}} = \gamma_{v_1}^{[2]}(1, \vec{v}_1^{[2]}) $$
$$ (v^{\mu}_2)^{\text{[2]}} = (1, \vec{0}) $$
(in the above and what follows, read superscript $[1]$ and $[2]$ to mean "measured in frame of particle 1/2").
Now simply compute $(1)$ in both frames, recalling that the results must be equal due to the invariance of said dot product. You will find quite easily that:
$$ \gamma_{v_1}^{[2]} = \gamma_{v_2}^{[1]}\\
\Rightarrow \frac{1}{\sqrt{1-\big|\big|\vec{v}_1^{[2]}\big|\big|^2}} = \frac{1}{\sqrt{1-\big|\big|\vec{v}_2^{[1]}\big|\big|^2}} $$
which clearly implies the magnitudes of $\vec{v}_1^{[2]}$ and $\vec{v}_2^{[1]}$ are equal. Now how can we also demonstrate that they are colinear but opposite in direction? To do that, consider what happens when you simply rotate the lab's coordinate system counterclockwise by an angle $\alpha$ such that:
$$ \sin\alpha = \frac{v_x}{\sqrt{v_x^2+v_y^2}}$$
$$ \cos\alpha = \frac{v_y}{\sqrt{v_x^2+v_y^2}}$$
I'll leave this for you as an exercise to demonstrate that in this rotated coordinate system (denoted by $(x',y')$ say), both particles move with a common identical velocity along the $x'$ axis and opposite (but generally not equal in magnitude!) velocities along the $y'$ axis. Therefore, if one further boosts into a frame that is itself translating along this rotated $x'$ axis with the appropriate velocity (details left as an exercise), the problem becomes one dimensional, and we can argue by symmetry that the rest frames of both particles must measure the other's velocity to be equal but opposite along that one dimensional line.
Finally, having reached this point, it should not be difficult to show that any additional rotation/boost of the coordinate system from this last one cannot introduce an asymmetry of the relative velocity as measured in one frame relative to the other (so that the conclusion is true even after we "unwind" the boost and rotation we've done).
To upshot of all that is that the entire structure of relativity, SR in particular, relies on the fact that the relative velocity of two inertial frames is a property of the transformations both from and to each of them. Therefore without this symmetry, or in other words the $\gamma$ factor being equal both in the "1 to 2" and "2 to 1" Lorentz transformations, SR wouldn't really make sense.
Now note that $\gamma = (1-\vec{v}^2)^{-1/2}$, so $\gamma$ can and will indeed be equal in both said rest frames of the particles, while $\vec{v}$ as measured in both will be equal in magnitude and opposite in direction.
Edit
I'll also try to address the following part of your question, because it appears to entail some conceptual loose ends where you may benefit from a bit more assistance:
Solving for velocity of 2 with respect to 1, let $V_2 = (0,v_y,0), V_1 = (v,0,0)$
$$V_{21} = (\frac{0-v}{1-\frac{0*v}{c^2}},\frac{v_y* \sqrt{1-\frac{v^2}{c^2}}}{1-\frac{0*v}{c^2}},0)$$
Note that as far as I can see, what you got is correct. You may be baffled by the asymmetry you get if you do the same to find $V_{12}$, to use a notation a bit more similar to yours, this will give you:
$$V_{12} = \left(\frac{v_x}{\gamma_{v_y}}, -v_y\right) $$
and then $V_{21}$ and $V_{12}$ as vectors are clearly not colinear. But there's no reason to expect them to be! Because you have not defined here a transformation directly from $1$ to $2$ or vice versa. Instead, when you computed the vector $V_{21}$ you essentially boosted from the lab frame to the frame of $1$ and found the three velocity of $2$. In what I've added, the result similarly comes from boosting from the lab frame to that of $2$ and finding the three velocity of $1$.
Since we have only thereby (implicitly) used two separate transformations, connecting either $1$ or $2$ to the lab coordinates, there's no reason that the relative velocity of $1$ found from the "Lab to 2" transformation will be equal and opposite as a vector to the velocity of $2$ found from the "Lab to 1" transformation. I hope that's clear and not just a word salad (if it is, read on, we'll try to make it more precise shortly).
But on the other hand, the magnitudes of $V_{21}$ and $V_{12}$ are still equal (hence this quantity is symmetric in $V_1$ and $V_2$) because as we argued before, by also rotating the coordinate system corresponding to the rest frame of each object, we can always make the relative velocity of the other lie along a single axis. Therefore there's no (sane) way in which the magnitude of the relative velocity can differ.
You can also convince yourself of this more explicitly. After boosting to the frame of $1$, construct the explicit Lorentz transformation from this frame to $2$, so first boost to rest frame of $1$:
$$ t' = \gamma_{v_x}(t_{\text{lab}}-v_x x_{\text{lab}}) $$
$$ x' = \gamma_{v_x}(x_{\text{lab}}-v_x t_{\text{lab}}) $$
$$ y' = y $$
then use the found velocity of $2$ wrt $1$ to construct a direct transformation:
\begin{align*}
t'' &= \gamma_{\small{V_{21}}}\ (t'+v_x x' - v_y \sqrt{1-v_x^2} y') \\
x'' &= \gamma_{\small{V_{21}}} v_x t' + \left(1+\frac{\gamma_{\small{V_{21}}}^2}{1+\gamma_{\small{V_{21}}}}v_x^2\right)x'-\left(\frac{\gamma_{\small{V_{21}}}^2}{1+\gamma_{\small{V_{21}}}}v_x v_y\sqrt{1-v_x^2}\right)y'\\
y'' &= -\gamma_{\small{V_{21}}} v_y\sqrt{1-v_x^2} t' -\left( \frac{\gamma_{\small{V_{21}}}^2}{1+\gamma_{\small{V_{21}}}}v_x v_y\sqrt{1-v_x^2}\right) x'+ \left(1+\frac{\gamma_{\small{V_{21}}}^2}{1+\gamma_{\small{V_{21}}}}v_y^2(1-v_x^2)\right)y'
\end{align*}
Ok, this may seem rather painful, because a Lorentz boost in a general direction kind of just is but it's very simple to check this gives what we want. To find the three velocity of $1$ now in the rest frame of $2$, just put $x'=y'=0$! This gives:
$$ t'' = \gamma_{\small{V_{21}}} t' $$
$$ x'' = \gamma_{\small{V_{21}}} v_x t' = v_x t''$$
$$ y'' = -\gamma_{\small{V_{21}}} v_y\sqrt{1-v_x^2}t' = -v_y\sqrt{1-v_x^2} t'',$$
where in the last two equations we have plugged in the first.
It's easy now to read off that in this frame:
$$ (\dot{x}'',\dot{y}'')_{12} = (v_x, -v_y\sqrt{1-v_x^2}) $$
equal and opposite to your $V_{21}$, as expected!
