A DNA polymer with hundreds or thousands base pairs is modeled as a thin elastic rod. To find the equilibrium configurations and associated elastic energies as a function of linking number difference (ΔLk), a finite element scheme based on Kirchhoff's rod theory is newly formulated so as to be able to treat self-contact. The numerical results obtained here agree well with those already published, both analytical and numerical, but a much more detailed picture emerges of the several equilibrium states which can exist for a given ΔLk. Of particular interest is the discovery of interwound states having odd integral values of the writhing number and very small twist energy. The existence of non-circular cross-section of the rod, inhomogeneous elastic constants, and intrinsic curvature can all be readily incorporated into the new formalism.