Simple circuit, mesh method, can't resolve
Your schematic and the current directions I chose to take in my own head:
simulate this circuit – Schematic created using CircuitLab
I get the same equations you do. My work is on the left, and yours (or closer to yours) is on the right:
$$\small\begin{align*} -I_1\:R_1 - \left(I_1-I_3\right)R_2+V_1-V_2 &=0\:\text{V} & I_1\:R_1 + \left(I_1-I_3\right)R_2 &=V_1-V_2\\\\ V_2-\left(I_2-I_3\right)R_4-I_2\left(R_5+R_6\right) &= 0\:\text{V} &\left(I_2-I_3\right)R_4+I_2\left(R_5+R_6\right)&=V_2\\\\ -V_1-\left(I_3-I_1\right)R_2-I_3\:R_3 -\left(I_3-I_2\right)R_4&=0\:\text{V} & \left(I_3-I_1\right)R_2+I_3\:R_3 +\left(I_3-I_2\right)R_4&=-V_1 \end{align*}$$
These result in the same matrix you formed:
$$\left[\begin{array}{ccc} R_1+R_2 & 0 & -R_2\\\\ 0 & R_4+R_5+R_6 & -R_4\\\\ -R_2 & -R_4 & R_2+R_3+R_4 \end{array}\right]\left[\begin{array}{c}I_1\\\\I_2\\\\I_3\end{array}\right]=\left[\begin{array}{c}V_1-V_2\\\\V_2\\\\-V_1\end{array}\right]$$
If you want to solve this by hand, Cramer's rule is often taught. Here:
$$\tiny\begin{align*} I_1 &= \frac{\text{det}\left[\begin{array}{ccc}V_1-V_2 & 0 & -R_2\\ V_2 & R_4+R_5+R_6 & -R_4\\ -V_1 & -R_4 & R_2+R_3+R_4\end{array}\right]}{\text{det}\left[\begin{array}{ccc}R_1+R_2 & 0 & -R_2\\ 0 & R_4+R_5+R_6 & -R_4\\ -R_2 & -R_4 & R_2+R_3+R_4\end{array}\right]}& I_2 &= \frac{\text{det}\left[\begin{array}{ccc}R_1+R_2 & V_1-V_2 & -R_2\\ 0 & V_2 & -R_4\\ -R_2 & -V_1 & R_2+R_3+R_4\end{array}\right]}{\text{det}\left[\begin{array}{ccc}R_1+R_2 & 0 & -R_2\\ 0 & R_4+R_5+R_6 & -R_4\\ -R_2 & -R_4 & R_2+R_3+R_4\end{array}\right]}& I_3 &= \frac{\text{det}\left[\begin{array}{ccc}R_1+R_2 & 0 & V_1-V_2\\ 0 & R_4+R_5+R_6 & V_2\\ -R_2 & -R_4 & -V_1\end{array}\right]}{\text{det}\left[\begin{array}{ccc}R_1+R_2 & 0 & -R_2\\ 0 & R_4+R_5+R_6 & -R_4\\ -R_2 & -R_4 & R_2+R_3+R_4\end{array}\right]} \end{align*}$$
So have you tried this approach?
A recent paper, A condensation-based application of Cramer’s rule for solving large-scale linear systems, Habgood & Arel, Journal of Discrete Algorithms 10, 2012, pp. 98–109, demonstrates that Cramer's rule can perform on the same computational order as other methods, such as LU decomposition. Nice to read about that.
Assuming you've worked out the three currents above, you can now easily work out the magnitude of the voltage across \$R_5\$ (if I understood your use of \$V_x\$, correctly) and you can work out its polarity with respect to the presumed current, \$I_2\$, as well.
I think you know how to calculate the current in \$V_1\$ as the sum of two of the three currents you worked out above. And similarly, for \$V_2\$. From those currents and the known voltages of each, you should have no trouble working out the power they contribute (the sign telling you if they are generating (-) or dissipating (+) power.)
These last details shouldn't be difficult, once you have the three currents worked out.
If you want to try here is the circuit:
Firstly, try and look at things as an EE would. For instance I could convert your diagram like this: -
Note that mentally I've named a node as 0 volts because it makes sense to me but it isn't an important thing to do if you don't want to do this. Then I've split V2 into two occurances (V2 and V3) so that I can disconnect the net between R1 and R6. I've also combined R5 and R6 into R7.
What I'd do next is make V1 and R2 into a 2 amp current source in parallel with R2 and make V3 and R1 into another 2 amp current source in parallel with R1.
Those two parallel current sources make 4 amps across 3.333 ohms (R1 || R2). I can then convert back to a common voltage source of 13.333 volts in series with 3.333 ohms and the job is looking a lot simpler: -
Next I would convert V2 and R7 into a current source of 1.818 amps in parallel with R7 and this then is placed in parallel with R4. Combine R4 with R7 and convert back to a voltage source and you are at the important point of calculating the current through R3. From hereon in it becomes a case of using R3's current to help solve the more localized voltages on nodes.