Why does MOM fail for large number of subdivision for finite dipole?

In Constantine A. Balanis' book about antennas, he introduced the method of moments for current distribution over a finite dipole.

I found that the method of moments works very bad for a half-wavelength dipole when you increase the subsegment to 2001. What's wrong?

POCKLINGTON''S EQN NUMBER OF SUBDIVISIONS <ODD NUMBER> = 2001 TOTAL DIPOLE LENGTH <WAVELENGTHS> = 0.9 RADIUS OF DIPOLE <WAVELENGTHS> =0.005 MAGNETIC-FRILL 

The following is what it looks like when NUMBER OF SUBDIVISIONS = 2001

The following is what it looks like when NUMBER OF SUBDIVISIONS = 201

The link for this code in matlab. get code here

1. Different Ways of Discretizing the Integral

Original Simpson’s-rule–style code

1. Splits the integral over $z'$ into segments and uses Simpson’s rule to approximate
   $$ \int_{-l/2}^{+l/2} \frac{e^{-j k R}}{R^5} \bigl[(1 + j k R)(2R^2 - 3a^2) + (k a R)^2\bigr] \,dz' $$

$$\int_{-l/2}^{+l/2} \frac{I_z(z') e^{-jkR}}{4\pi R^5} \left[ (1 + jkR)(2R^2 - 3a^2) + (kaR)^2 \right] dz' = -j\omega E_z^i (\rho = a) $$ 2. There are explicit factors like $\tfrac{4}{3}$, or sums of coefficients 4 and 2 in the loop (if mod(k,2)~=0 ...), which is the Simpson integration logic.

Matrix-based Method of Moments (MoM) code

1. Builds a discrete matrix $\mathbf{Z}$ (which you named imp) of dimension $N\times N$:
   $$ Z_{nm} \;=\; f\bigl(a,z_n, z_m\bigr)\,\Delta z' \quad\text{where}\quad z_n,\;z_m \;\text{are discrete points along the wire, etc.} $$
2. The function

   function zi = f(a,z,z1) R = sqrt( (z - z1)^2 + a^2 ); k = 2*pi; p1 = exp(-1i*k*R)/(4*pi*R^5); p2 = (1+1i*k*R)*(2*R^2 - 3*a^2) + (k*a*R)^2; zi = p1*p2; end 

   corresponds to the kernel of your integral $\frac{1}{4\pi R^5}\bigl(\cdots\bigr)$.
3. Instead of using Simpson’s rule, the code essentially multiplies by $\Delta z'$ (dz1) and then sums these entries in the matrix equation $\mathbf{Z}\,\mathbf{I} = \mathbf{V}$.

Key difference: The Simpson’s-rule code accumulates the integral by adding up function values with coefficients (4,2,1). The matrix code lumps the kernel values into matrix entries (Z_{n,m}) (plus a uniform $\Delta z'$ factor), then solves $\mathbf{Z}\,\mathbf{I} = \mathbf{b}$.

Constantine A. Balanis used simpson rule in his matlab code.