What is the closed-loop gain in this op amp circuit?

Assuming this is a problem from a book and we can assume ideal parameters to the opamp where no specifications are otherwise given, then:

^{simulate this circuit – Schematic created using CircuitLab}

Assuming the equivalent circuit at the top for the opamp, the nodal equations are:

$$\begin{align*} \frac{V_X}{R_1}+\frac{V_X}{R_2}+\frac{V_X}{R_{IN}}&=\frac{V_O}{R_1}+\frac{V_S}{R_{IN}}\\\\ \frac{V_O}{R_1}+\frac{V_O}{R_3}+\frac{V_O}{R_{OUT}}&=\frac{V_X}{R_1}+\frac{\left(V_S-V_X\right)\cdot A_\text{OL}}{R_{OUT}}\\\\\therefore\\\\V_X&\approx 0.999955\cdot V_S\\\\V_O&\approx 8.999593\cdot V_S \end{align*}$$

The script I used in sympy (worth getting) is:

var('r1 r2 r3 ri ro aol vo vx vs')
e1=Eq(vx/r1+vx/r2+vx/ri,vo/r1+vs/ri)
e2=Eq(vo/r1+vo/r3+vo/ro,vx/r1+(vs-vx)*aol/ro)
ea=solve([e1,e2],[vx,vo])
    { vo: r1*r3*vs*(aol*(r1*r2 + r1*ri + r2*ri) - r2*(aol*r1 - ro))/(r2*r3*ri*(aol*r1 - ro) + (r1*r2 + r1*ri + r2*ri)*(r1*r3 + r1*ro + r3*ro)),
      vx: r1*r2*vs*(aol*r3*ri + r1*r3 + r1*ro + r3*ro)/(r2*r3*ri*(aol*r1 - ro) + (r1*r2 + r1*ri + r2*ri)*(r1*r3 + r1*ro + r3*ro))
    }
ea[vx].subs({aol:2e5,ri:2e6,ro:50,r1:40e3,r2:5e3,r3:20e3})
0.999954839544092*vs
ea[vo].subs({aol:2e5,ri:2e6,ro:50,r1:40e3,r2:5e3,r3:20e3})
8.99959265268771*vs

Just nodal is enough here. The main thing is to figure out the opamp model from the specs you were given.

In this type of circuit, one easy way is to apply superposition and determine the control variable, \$\epsilon\$. To express it, we will apply superposition to this circuit featuring 1 controlled source: first we consider \$V_{in}=0\$ and second, we will consider the op-amp output \$\epsilon A_{OL}\$ equal to 0. In the first case, the circuit is below:

If you do the maths ok, you can determine \$\epsilon_1\$ whose value is in the Mathcad sheet at the end. Now, put \$V_{in}\$ back in place and consider \$\epsilon A_{OL}\$ equals 0:

Again, if you do the maths ok, you can determine \$\epsilon_2\$ whose value is in the Mathcad sheet at the end. You have the final \$\epsilon\$ value by writing \$\epsilon=\epsilon_1+\epsilon_2\$ and solve for \$\epsilon\$. That is what Mathcad did for us.

Now that you are there, you are almost done. Looking at the circuit, you see that \$V_{out}=\epsilon A_{OL} - i_1R_o\$. \$i_1\$ is the sum of two currents that you easily determine. You end-up with an equation featuring \$V_{out}\$, \$V_{in}\$ and \$\epsilon\$ that you now have on hand. If Mathcad does the job ok, for a 1-V input voltage, the output is exactly 8.99183 V for an open-loop gain of 10k:

The Mathcad file is here but I won't simplify equations, it's already quite late over here : )

These techniques using superposition are part of the fast analytical circuits techniques or FACTs. A solving technique that I encourage students and EEs to acquire.

Edit: I did not realize that I did not use the exact values given in the post. If I plug these in the sheet, I have a gain of 8.99959 as Monsieur jonk found.