I used MATLAB eig() to find eigenvectors and eigenvalues of a complex symmetric matrix. I searched through MATLAB online documentation to find a link to some algorithm, but failed. My curiosity is piqued also because of the fact that the algorithm used by eig() didn't seem to be something simple enough. I am saying this because we have a rudimentary conjugate gradient complex symmetric eigensolver in FORTRAN, and we get poor quality of complex orthogonality* between eigenvectors, unlike MATLAB.
*note that for a complex symmetric matrix, eigenvectors corresponding to distinct eigenvalues have a zero transpose inner product, not a zero conjugate-transpose inner product. That is, $v_1^Tv_2=0$, but $v_1^\dagger v_2\neq 0$.
The exact algorithm isn't known because Matlab tweaks it however it is based on LAPACK which can be found online at the Netlib. They are available in Fortran. If you read it is variants of the shifted QR algorithm.
The routines are here.