To study linearized magnetohydrodynamic (MHD) waves, continuous spectra, and instabilities in coronal magnetic flux tubes that are anchored in dense chromospheric and photospheric regions, a two-dimensional numerical code, called PARIS, has been developed. PARIS solves the pertinent nonlinear Grad-Shafranov type, partial differential equation for the magnetic flux on a flux coordinate grid. Both a straight field line coordinate system and an orthogonal flux coordinate system are exploited. Isoparametric bicubic Hermite finite elements have been adopted to solve the Grad-Shafranov-like equation. These elements allow for a continuous representation of the Bur and the gradient of the flux throughout the tube and can be aligned conveniently along the boundary of the tube. These properties are important to obtain an accurate representation of the solution on flux coordinate grids. An analytical test case is used to show that accurate solutions have been obtained, even for a small number of grid points. The equilibria calculated by PARIS are used to study the continuous spectra of two-dimensional magnetic flux tubes. One illustrative example is given here; extensive results are presented elsewhere (A.J.C. Belien, S. Poedts and J.P. Goedbloed, Astron. Astrophys. 322 (1997) 995). The equilibria obtained by PARIS are also well suited to calculate the stability and the normal mode MHD spectrum. (C) 1997 Elsevier Science B.V.