It's probably better to use the dependent voltage source, as you did. But I'd infer it mentally and keep the schematic more like this:
simulate this circuit – Schematic created using CircuitLab
This allows me to quickly set up four equations and four unknowns. The purpose of \$I_{_\text{O}}\$ is to allow me to inject a current. I'd set it to \$0\:\text{A}\$ and then \$1\:\text{A}\$ and measure the difference of \$V_{_\text{O}}\$ (and divide by 1, the obvious change in current.)
Keeping all this in symbolic form can be handled easily with SymPy. So that's my recommendation.
Please note that in your paper-written diagram, you've not labeled the node between \$R_1\$ and \$R_2\$. I don't know if that's an oversight or if you are fully aware, but just didn't say anything about it. And I don't see how you can move towards the solution you show without it clearly named.
var( 'vi vo va ia r1 r2 rout vm avo io' ) # list needed variables
eq1 = Eq( vm/r1 + vm/r2, vi/r1 + vo/r2 ) # KCL for VM
eq2 = Eq( va, -avo*vm ) # opamp open loop voltage gain
eq3 = Eq( vo/rout + vo/r2, va/rout + v/r2 + io ) # KCL for VO, with injection IO current
eq4 = Eq( va/rout, vo/rout + ia ) # KCL for VA, with opamp output current
ans = solve( [ eq1, eq2, eq3, eq4 ], [ v, ia, vo, va ] )
vo0 = ans[vo].subs( { io:0 } ) # VO without injected current
vo1 = ans[vo].subs( { io:1 } ) # VO with injected 1A current
pprint( simplify( vo1 - vo0 ) ) # Print output resistance
rout⋅(r₁ + r₂)
───────────────────────
avo⋅r₁ + r₁ + r₂ + rout
The above is exactly equivalent to your given correct answer. This can be tested in the following way:
n = rout*(r1+r2)/(r1+r2+rout) # ROUT || (R1 + R2)
d = 1 + r1*avo/(rout+r1+r2) # 1 + R1*Avo/(ROUT + R1 + R2)
pprint( simplify( n / d ) )
rout⋅(r₁ + r₂)
───────────────────────
avo⋅r₁ + r₁ + r₂ + rout
Their answer is correct and the approach shown above is soundly reasoned in getting to the same place.