I don't see how you get Iin = VxGd from the 2nd equation. I agree with you that the two node voltages are proportional from the 1st eqn, but when I solve for Iin in the 2nd eqn (eliminating Vx via your first equation) I get a term -gds(gds-g)/(gds+GD), a partial current, that must be flowing through the independent Vgs source. Applying KCL at the botttom node, then seems to contradict your later statement that Iin = Iout.
The transfer function is: $$\small G(s)=\frac{1}{(RC)^2s^2+3RCs+1}$$
Hence \$\omega_n=\frac{1}{RC}=\small10^4\:rad/s\:(=1592\:Hz)\$, and \$\small\zeta=1.5\$, and it can be seen that the DC gain (\$\small s=0\$) is unity.
Converting this to the frequency domain, using \$ s\rightarrow j\omega\$:
$$\small G(j\omega)=\frac{1}{1-(\omega RC)^2+j3\omega RC}$$
At \$\small 10\:\small Hz\$, \$\small \omega RC=0.00628\$, hence the gain is almost unity and the phase angle is almost zero. At \$\small 1\:\small kHz\$, \$\small \omega RC=0.628\$, giving a gain of \$\small 0.505\$, and phase angle of \$\small \phi=-72^o\$.
So it seems that there's a problem with your experimental set-up. What's the input impedance of the instrument measuring Vout?
Let's do some detective work:
If the input impedance of the instrument were \$\small 3 \: k\Omega\$ resistive, then (i) the gain and phase at DC would be \$\small 0.6\$ and zero, respectively (i.e. same as your results); and (ii) the gain and phase at \$\small 1590 \:\small Hz\$ would be \$\small 0.29\$ and \$\small -79^o\$, which compares with your measurement of \$\small 0.31\$ and \$\small -73^o\$.
Best Answer
I think that writing the loop equations would be easier.
The loop equations for first two loops: $$I_1(Z+R) = V_{in}$$ $$- I_1R +I_2(Z+2R) -I_3R = 0$$
Where \$Z=\dfrac{1}{Cs}\$. From this:
$$I_2(Z+2R) -I_3R = V_{in}\frac{R}{R+Z}\tag1$$ The remaining two loop equations: $$- I_2R +I_3(Z+2R) -I_4R = 0\tag2$$ $$- I_3R +I_4(Z+2R) = 0\tag3$$
Expressing in matrix form:
$$\left[\begin{array}{ccc} &Z+2R &-R &0 \\ &-R &Z+2R &-R \\ &0 &-R & Z+2R \end{array}\right] \left[\begin{array}{c} I_2\\ I_3\\ I_4 \end{array}\right] = \left[\begin{array}{c} \frac{V_{in}R}{R+Z}\\ 0\\ 0 \end{array}\right]$$
Now by Cramer's rule: $$ I_4 = \frac{\left|\begin{array}{ccc} Z+2R &-R &\frac{V_{in}R}{R+Z} \\ -R &Z+2R &0 \\ 0 &-R & 0 \end{array}\right|}{\left|\begin{array}{ccc} Z+2R &-R &0 \\ -R &Z+2R &-R \\ 0 &-R & Z+2R \end{array}\right|}$$
$$V_{out} = I_4\times R$$
From this the transfer function can be calculated. Gain and phase shift can be calculated from transfer function. (substitute \$Z=\frac{1}{jwC}\$)