Assume your \$V_1\$ voltage source is called \$V_X\$ and that the current in \$V_X\$ is \$I_X\$, with the arrow pointing from \$x_9\$ to its tip at \$x_8\$. Then:
$$\begin{align*}
\frac{V_1}{R_1}+\frac{V_1}{R_3}&=\frac{V_2}{R_1}+\frac{V_4}{R_3}\tag{x1}\\\\
\frac{V_2}{R_1}+\frac{V_2}{R_5}+\frac{V_2}{R_9}&=\frac{V_1}{R_1}+\frac{V_5}{R_5}+\frac{V_3}{R_9}\tag{x2}\\\\
\frac{V_3}{R_9}+\frac{V_3}{R_{10}}&=\frac{V_2}{R_9}+\frac{V_6}{R_{10}}\tag{x3}\\\\
\frac{V_4}{R_3}+\frac{V_4}{R_4}+\frac{V_4}{R_6}&=\frac{V_1}{R_3}+\frac{V_5}{R_4}+\frac{V_7}{R_6}\tag{x4}\\\\
\frac{V_5}{R_4}+\frac{V_5}{R_5}+\frac{V_5}{R_8}+\frac{V_5}{R_{12}}&=\frac{V_4}{R_4}+\frac{V_2}{R_5}+\frac{V_8}{R_8}+\frac{V_6}{R_{12}}\tag{x5}\\\\
\frac{V_6}{R_{10}}+\frac{V_6}{R_{11}}+\frac{V_6}{R_{12}}&=\frac{V_3}{R_{10}}+\frac{V_9}{R_{11}}+\frac{V_5}{R_{12}}\tag{x6}\\\\
\frac{V_7}{R_6}+\frac{V_7}{R_7}&=\frac{V_4}{R_6}+\frac{V_8}{R_7}\tag{x7}\\\\
\frac{V_8}{R_7}+\frac{V_8}{R_8}&=I_X+\frac{V_7}{R_7}+\frac{V_5}{R_8}\tag{x8}\\\\
V_9&=0\:\text{V}\tag{x9}\\\\
V_8&=V_9+V_X\tag{10}\\\\
\end{align*}$$
The above is 10 equations and 10 unknowns. You should be able to populate the matrices from the above. To do that, just move terms around so that your variables are all on the left side with associated factors and there are only constants on the right side. (The factors on the left will then be elements in your left-side matrix.)
There are some obvious simplifications that would allow you to reduce the unknowns and the equations. But I'm sticking with your approach, labeling all those nodes as you did regardless that you could have simplified things a bit.
Since someone has seen fit to mark this answer wrong (-1), I'll take the time to prove that it produces accurate and reliable results. Too bad they were too cowardly to identify themselves and their problem. But I can at least show them their error in thinking about what I wrote.
In Sage, this is how you might formulate the above:
var('v1 v2 v3 v4 v5 v6 v7 v8 v9 ix vx')
var('r1 r2 r3 r4 r5 r6 r7 r8 r9 r10 r11 r12')
eq1= Eq( v1 * (1/r1 + 1/r3), v2/r1 + v4/r3 )
eq2= Eq( v2 * (1/r1 + 1/r5 + 1/r9), v1/r1 + v5/r5 + v3/r9 )
eq3= Eq( v3 * (1/r9 + 1/r10), v2/r9 + v9/r10 )
eq4= Eq( v4 * (1/r3 + 1/r4 + 1/r6), v1/r3 + v5/r4 + v7/r6 )
eq5= Eq( v5 * (1/r4 + 1/r5 + 1/r8 + 1/r12), v4/r4 + v2/r5 + v8/r8 + v6/r12 )
eq6= Eq( v6 * (1/r10 + 1/r11 + 1/r12), v3/r10 + v9/r11 + v5/r12 )
eq7= Eq( v7 * (1/r6 + 1/r7), v4/r6 + v8/r7 )
eq8= Eq( v8 * (1/r7 + 1/r8), ix + v7/r7 + v5/r8 )
eq9= Eq( v9 , 0 )
eq10=Eq( v8 , v9 + vx )
ans=solve([eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,eq10],[v1,v2,v3,v4,v5,v6,v7,v8,v9,ix])
[(i, ans[i].subs({r1:.5,r3:.1,r4:.1,r5:.5,r6:.2,r7:.2,r8:.1,r9:.2,r10:.1,r11:.2,r12:.5,vx:9})) for i in ans]
[(v9, 0),
(v4, 7.67171239356670),
(ix, 16.4545884578997),
(v6, 3.29091769157994),
(v3, 4.05723746452223),
(v8, 9),
(v7, 8.33585619678335),
(v2, 5.58987701040681),
(v1, 7.32473982970672),
(v5, 7.68661305581835)]
The above uses a simple loop to walk through the list and emit the values, but in no particular order.
If you want them sorted, you could instead write:
[(eval(key), ans2[eval(key)].subs({r1:.5,r3:.1,r4:.1,r5:.5,r6:.2,r7:.2,r8:.1,r9:.2,r10:.1,r11:.2,r12:.5,vx:9})) for key in sorted([str(key) for key in ans2])]
[(ix, 16.4545884578997),
(v1, 7.32473982970672),
(v2, 5.58987701040681),
(v3, 4.05723746452223),
(v4, 7.67171239356670),
(v5, 7.68661305581835),
(v6, 3.29091769157994),
(v7, 8.33585619678335),
(v8, 9),
(v9, 0)]
You can also hand-request them like this:
ans[v1].subs({r1:.5,r3:.1,r4:.1,r5:.5,r6:.2,r7:.2,r8:.1,r9:.2,r10:.1,r11:.2,r12:.5,vx:9})
7.32473982970672
ans[v2].subs({r1:.5,r3:.1,r4:.1,r5:.5,r6:.2,r7:.2,r8:.1,r9:.2,r10:.1,r11:.2,r12:.5,vx:9})
5.58987701040681
ans[v3].subs({r1:.5,r3:.1,r4:.1,r5:.5,r6:.2,r7:.2,r8:.1,r9:.2,r10:.1,r11:.2,r12:.5,vx:9})
4.05723746452223
ans[v4].subs({r1:.5,r3:.1,r4:.1,r5:.5,r6:.2,r7:.2,r8:.1,r9:.2,r10:.1,r11:.2,r12:.5,vx:9})
7.67171239356670
ans[v5].subs({r1:.5,r3:.1,r4:.1,r5:.5,r6:.2,r7:.2,r8:.1,r9:.2,r10:.1,r11:.2,r12:.5,vx:9})
7.68661305581835
ans[v6].subs({r1:.5,r3:.1,r4:.1,r5:.5,r6:.2,r7:.2,r8:.1,r9:.2,r10:.1,r11:.2,r12:.5,vx:9})
3.29091769157994
ans[v7].subs({r1:.5,r3:.1,r4:.1,r5:.5,r6:.2,r7:.2,r8:.1,r9:.2,r10:.1,r11:.2,r12:.5,vx:9})
8.33585619678335
ans[v8].subs({r1:.5,r3:.1,r4:.1,r5:.5,r6:.2,r7:.2,r8:.1,r9:.2,r10:.1,r11:.2,r12:.5,vx:9})
9
ans[v9].subs({r1:.5,r3:.1,r4:.1,r5:.5,r6:.2,r7:.2,r8:.1,r9:.2,r10:.1,r11:.2,r12:.5,vx:9})
0
ans[ix].subs({r1:.5,r3:.1,r4:.1,r5:.5,r6:.2,r7:.2,r8:.1,r9:.2,r10:.1,r11:.2,r12:.5,vx:9})
16.4545884578997
You will note that these values (printed using any method) will be precisely matched by any Spice simulation you composite up. As I wrote earlier, this approach does work. Though of course there are simplified approaches that would get the more important results from which you could easily get the rest.
My personal opinion for folks caring about symbolic matrix solutions to electronic problems is to become very familiar with SageMath and Python. These are incredibly powerful and free tools available to anyone and well worth a small investment into their learning curves.