How to find the Q point of the NMOS transistor in Voltage-Divider Biasing Circuit?
As you've noted

$$V_{GS} = V_G - I_D R_S = 3.33 - (3000\Omega)I_D$$

Rearranging:

$$I_D = \frac{3.33\text{V} - V_{GS}}{3000\Omega} \tag1$$

Your two equations for \$I_D\$ correspond to the MOSFET in the linear and saturation regions, respectively. Only one of the two applies to this circuit, depending on whether the MOSFET is in the linear or saturation regions. Guess that the MOSFET is in saturation, which means

$$I_D = K_n(V_{GS} - V_{TN})^2 \tag2$$

Set \$(1)\$ and \$(2)\$ equal to each other and solve for \$V_{GS}\$ (the only unknown):

$$\frac{3.33\text{V} - V_{GS}}{3000\Omega} = K_n(V_{GS} - V_{TN})^2$$

You will get two solutions: \$V_{GS} = -.623\$V or \$V_{GS} = 1.9568\$V. The first solution makes no sense, so \$V_{GS} = 1.9568\$V. Plug this back into \$(1)\$ or \$(2)\$ to find \$I_D = 458\mu\$A.

Now we need to check that the MOSFET is really in saturation so that \$(2)\$ is the right equation to use (rather than your first equation for \$I_D\$). In saturation

$$V_{DS} &gt; V_{GS} - V_{TN}$$

We have \$V_D = 10\text{V} - I_D R_D = 8.17\$V and \$V_S = I_D R_S = 1.37\$V so

$$V_{DS} = V_D - V_S = 6.8\text{V} &gt; V_{GS} - V_{TN} = 1.9568\text{V} - 1\text{V} = 0.9578\text{V}$$

so the MOSFET is indeed in saturation and our solutions for \$I_D\$, \$V_{GS}\$ and \$V_{DS}\$ are correct.

If you simulate your schematic in CircuitLab (I've edited it so you can simulate it with the given MOSFET \$V_{TN}\$ and \$K_n\$) you'll find that it gives \$I_D \approx 380\mu\$A, which is not far from my solution (the slight discrepancy is due to other MOSFET simulation parameters, which are unknown).