"Which brings me to the main point I'm missing: How to find out where those added poles are? For example, how did you calculated the added pole due to Cgd?"
Assuming Cgd is 500pF (which can be seen by looking at the IRF3205's datasheet; look in Fig.5 at the plot for Crss—this is Cgd, the MOSFET's internal capacitance from gate to drain), the added pole occurs at 1/(2*pi*330Ω*500pF)=1MHz. Remembering the series RC we considered earlier, but replacing R=330Ω and C=10nF, we now place a capacitor in parallel to these of value 500pF. The impedance of this combo is very high at low frequencies and decreases until it reaches 50kHz, at which point it levels off at 330Ω, and then above 1MHz it falls again.
You can visualize this easily using LTspice. Here's a graph of the impedance, and the circuit that produced it:
I'm feeding the network a current of amplitude 1, so the voltage across the network represents its impedance (V=IZ; for I=1, V=Z). As you can see, with Cgd absent the impedance remains resistive at high frequencies (red), which was the intended operation. However, with the large Cgd introduced, the impedance becomes capacitive again above ~1MHz (orange), which brings the phase down and causes instability if the unity gain frequency occurs above it.
Looking at the orange graph, we notice three distinct regions. Starting from the left the impedance is falling, indicating a capacitive impedance (since for a capacitor, Z=1/sC). The fall is at 20dB/dec, which also indicates a single-pole rolloff (the pole occurred at zero frequency). Whenever you have a pole, it is accompanied by a falling gain—20dB/dec per pole. At 50kHz a zero is introduced by the resistor. If there weren't any poles, the zero would actually cause a 20dB/dec rise, but since there was already a pole it simply cancels out the fall, resulting in a level region above 50kHz. At 1MHz a second pole occurs, causing the impedance to fall again.
Doing the math explicitly, we have a parallel combination of (1/sCgd) // (R+1/sC). This translates to:
\[Z=\frac{1}{sCgd+sC/(sCR+1)}\].
\[=\frac{sCR+1}{sCgd(sCR+1)+sC}\].
\[=\frac{sCR+1}{s[sCgdCR+Cgd+C]}\].
\[=\frac{sCR+1}{\frac{s}{C+Cgd}[s\frac{CCgd}{C+Cgd}R+1]}\].
\[=\frac{sCR+1}{\frac{s}{C+Cgd}[s(C//Cgd)R+1]}\].
This equation expresses the impedance as a function of frequency, s (where s is in general a complex number, but we usually substitute jω for sinusoidal signals; I prefer to write s because it saves space). Occurrences of s in the numerator cause "zeros" in the frequency response (here, you can see that the impedance goes to zero if s=-1/CR), while occurrences in the denominator cause "poles" in the response (you can see the impedance goes to infinity if s=0 or if s=-1/(C//Cgd)R). For sinusoidal signals we won't actually see zeros and infinities here, but we will still see something distinctive happen—we'll see a change in the slope. We have a pole at ω=0 causing the impedance to fall from the start; then we have a zero at ω=1/CR causing the impedance to level off, and finally another pole at ω=1/(C//Cgd)R causing the impedance to fall again. Intuitively, if I'm graphing something in the numerator (such as y=x) then it should go up as I traverse the horizontal axis, and if I graph something in the denominator (such as y=1/x) it should go down. The more terms in the numerator (y=x² or x³), the faster it rises, and the more terms in the denominator (y=1/x² or 1/x³), the faster it falls. Same thing happens here.
One thing to notice is that for the second pole at ω=1/(C//Cgd)R, C»Cgd (10nF»500pF). In "parallel" combinations, the smaller value dominates, so we can approximate the pole frequency to ω=1/CgdR. This is in line with what was stated earlier, indicating that we didn't actually need to do all this math, although it's useful for building our confidence.
Stepping back for a bit, this all comes into play by affecting the MOSFET's gain. Its voltage gain in this configuration is (again, approximately) Zfeedback/Zdrive. Since Zdrive is a resistor of 6.8kΩ, this stage behaves as an integrator below 50kHz, transitions to an amplifier of gain 0.05, and then behaves like an integrator again above 1MHz. The gain curve will roughly follow the feedback impedance curve. The new R and C value don't guarantee stability when you add fifty more MOSFETs on, but it should have good margin with one.
As crutschow recommended, it's a good idea to simulate the circuit. You can find
a model of the IRF3205 from International Rectifier (hopefully its model is suitable; you can compare it to the datasheet for a double-check. Instead of reading the values, you can build simulations that try to mimic the datasheet graphs, and then see how close they are—this saves you the headache of trying to figure out what every single simulation parameter is). When simulating the circuit, this is what I get:
Obviously, I couldn't be bothered with figuring out how to get a nice symbol working so I added the MOSFET's subcircuit with a SPICE line. With any luck, this MOSFET amplifier shouldn't make the global loop unstable. Plug in a model for an op-amp, and evaluate the loop gain. This is what I get with an ideal, single-pole opamp with a 1MHz unity gain frequency (some call it the gain-bandwidth product, GBWP, or GBW):
This simulation shows about 45 degrees of phase margin at 59kHz, which most people consider acceptable. Going through the exercise of simulating sections and subsections, and tweaking values to see how things change, can really be enlightening. Simulations can also help catch issues you didn't realize would be there.