The exact within potential flow integral equation approach of Evans and Ford [Proc. R. Soc. London Ser. A 452, 373 (1996)] for the normal solitary wave, is here generalized to 2-layer, ‘‘internal’’ solitary waves. This differs in its mathematical form from other exact integral equation methods based on the complex velocity potential. For both ‘‘rigid lid’’ (i.e., flat toplayer surface) and ‘‘free-surface’’ boundary conditions, a set of coupled non-linear integral equations are derived by an application of Green’s theorem. For each point on the layer interface(s), these describe functional constraints on the profiles and interface fluid velocity moduli; the exact profiles and velocities being those forms that satisfy these constraints at all such interface points. Using suitable parametric representations of the profiles and interface velocity moduli as functions of horizontal distance, x, and utilizing tailored quadrature methods [Int. J. Comput. Math. B 6, 219 (1977)], numerical solutions were obtained by the Newton–Raphson method that are highly accurate even at large amplitudes. For ‘‘rigid lid’’ boundary conditions, internal wave solutions are presented for layer density and depth ratios typical of oceanic internal wave phenomena as found in the Earth’s marginal seas. Their various properties, i.e., mass, momentum, energy, circulation, phase and fluid velocities, streamline profiles, internal pressures, etc., are evaluated and compared, where possible, with observed properties of such phenomena as reported, for example, from the Andaman Sea. The nature of the limiting (or ‘‘maximum’’) internal wave is investigated asymptotically and argued to be consistent with two ‘‘surge’’ regions separating the outskirts flow from a wide mid-section region of uniform ‘‘conjugate flow’’ as advocated by Turner and Vanden-Broeck [Phys. Fluids 31, 286 (1988)].