Nodal Analysis algebra simplification
What I recommend doing is finding common denominator before plugging in, like this:

Node A

$$\frac{V_1-V_{ref}}{R_2} + \frac{V_1-V_{x}}{R_1} + \frac{V_1-V_{2}}{R_f} = 0 $$
$$ V_1(R_1R_2 + R_1R_f + R_2R_f) - V_2R_1R_2 - V_x(R_fR_2) + V_{ref}R_fR_1 = 0$$
$$V_x(R_fR_2) = V_1(R_1R_2 + R_1R_f + R_2R_f) - V_2R_1R_2 + V_{ref}R_fR_1$$

Node B

$$ \frac{V_2-V_{x}}{R_3} + \frac{V_2-V_{1}}{R_f} + \frac{V_2-V_{out}}{R_4} = 0 $$
$$ V_2(R_3R_4 + R_4R_f + R_3R_f) - V_1R_3R_4 - V_x(R_fR_4) + V_oR_fR_3 = 0$$
$$ V_x(R_fR_4) = V_2(R_3R_4 + R_4R_f + R_3R_f) - V_1R_3R_4 + V_oR_fR_3$$

since $$R_4=R_2 \space and \space R_3=R_1$$
$$=  V_x(R_fR_2) = V_2(R_1R_2 + R_2R_f + R_1R_f) - V_1R_1R_2 + V_oR_fR_1$$

Then plug it in:
$$V_1(R_1R_2 + R_1R_f + R_2R_f) - V_2R_1R_2 + V_{ref}R_fR_1 = V_2(R_1R_2 + R_2R_f + R_1R_f) - V_1R_1R_2 + V_oR_fR_1$$
$$V_oR_fR_1 = V_1(R_1R_2 + R_1R_f + R_2R_f) + V_1R_1R_2  - V_2(R_1R_2 + R_2R_f + R_1R_f) - V_2R_1R_2 + V_{ref}R_fR_1$$
I'm sure you can figure out the rest.