Been a while. But I finally see this. I'll give a hack at it. See if it sings.
Fast Attack
First off, redraw the schematic:

simulate this circuit – Schematic created using CircuitLab
This is so much simpler to analyze. All you need to do is Thevenize the left and right legs to get:

simulate this circuit
And from this you know that for maximum power then \$R\$ must be \$25\:\Omega\$. (Load resistance equals source resistance.)
It's actually pretty trivial to find.
Nodal
Let's go back to your schematic and do this the hard way. From what little I can glean from your writing, you have \$V=+40\:\text{V}\$ and you've assigned the voltage supply's negative terminal to \$0\:\text{V}\$. Let's look at what I think you started with:

simulate this circuit
Pretty much the same thing I wrote earlier. Just lots uglier.
To solve for the resistance of the known part of the circuit (the part seen by \$R\$ looking backwards into the circuit), you want the open-circuit voltage (\$V_{_\text{OC}}\$) that happens when \$R=\infty\:\Omega\$ and you want the short circuit current (\$I_{_\text{SC}}\$) through \$R\$ that happens when \$R=0\:\Omega\$. Then you can find that the maximum power value for \$R\$ happens when \$R=\frac{V_{_\text{OC}}}{I_{_\text{SC}}}\$.
With \$R=\infty\:\Omega\$ then:
eq1 = Eq( v/r1 + v/r3 + iv, v1/r1 + v2/r3 ) # KCL for node v
eq2 = Eq( v1/r1 + v1/r2, v/r1 + 0/r2 ) # KCL for node v1
eq3 = Eq( v2/r3 + v2/r4, v/r3 + 0/r4 ) # KCL for note v2
ans = solve( [ eq1, eq2, eq3 ], [ v1, v2, iv ] ) # Solve!!!
ans[v1].subs( { v:40, r1:80, r2:20, r3:10, r4:90 } )
8
ans[v2].subs( { v:40, r1:80, r2:20, r3:10, r4:90 } )
36
I'm using SymPy and Sage here. (All free to have.)
The difference is obviously \$\mid\, V_{_\text{OC}}\mid\:=28\:\text{V}\$. If we take the difference, v1-v2, then we find that \$V_{_\text{OC}}=v_1-v_2=-28\:\text{V}\$.
Now, let's compute the short-circuit voltages for both nodes. We'll keep the nodes distinct, but assign them equal to each other as that is what must be if we short-circuit the nodes. I've added a new current (ic) as shown below that proceeds from v1 to v2:
eq1 = Eq( v/r1 + v/r3 + iv, v1/r1 + v2/r3 ) # KCL for node v
eq2 = Eq( v1/r1 + v1/r2 + ic, v/r1 + 0/r2 ) # KCL for node v1
eq3 = Eq( v2/r3 + v2/r4, ic + v/r3 + 0/r4 ) # KCL for note v2
eq4 = Eq( v1, v2 ) # v1 = v2
ans = solve( [ eq1, eq2, eq3, eq4 ], [ v1, v2, iv, ic ] ) # Solve!!!
ans[ic].subs( { v:40, r1:80, r2:20, r3:10, r4:90 } )
-28/25
There we are.
So we need to compute \$R=\frac{V_{_\text{OC}}}{I_{_\text{SC}}}=\frac{-28\:\text{V}}{-\frac{28}{25}\:\text{A}}=25\:\Omega\$. That's the source resistance. And again, we know that \$R\$ must be the same. So \$R=25\:\Omega\$, as well, to maximize the power into \$R\$.
If you want to test all this, let's do a new set of nodal equations that now include \$R\$:
eq1 = Eq( v1/r1 + v1/r2 + v1/r, v/r1 + 0/r2 + v2/r ) # KCL for node v1
eq2 = Eq( v2/r3 + v2/r4 + v2/r, v/r3 + 0/r4 + v1/r ) # KCL for note v2
ans = solve( [ eq1, eq2 ], [ v1, v2] ) # Solve!!!
(ans[v1]-ans[v2]).subs( { v:40, r1:80, r2:20, r3:10, r4:90, r:25 } )
-14
From this, we'd expect to see the power as \$\frac{\left(-14\:\text{V}\right)^2}{25\:\Omega}=7.84\:\text{W}\$.
Let's see what LTspice says:

Yup! Looks good.
Feel free to slightly adjust \$R\$ up a little bit or down a little bit to see what the resulting power is. You will find that either way it gets a little smaller. The peak power remains only when \$R=25\:\Omega\$.
Nodal works. It just can be a pain at times.