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.