A computational study of superconducting states near the superconducting-normal phase boundary in mesoscopic finite cylinders is presented. The computational approach uses a finite element method to find numerical solutions of the linearized Ginzburg-Landau equation for samples with various sizes, aspect ratios, and cross-sectional shapes, i.e., squares, triangles, circles, pentagons, and four star shapes. The vector potential is determined using a finite element method with two penalty terms to enforce the gauge conditions that the vector potential is solenoidal and its normal component vanishes at the surface(s) of the sample. The eigenvalue problem for the linearized Ginzburg-Landau equations with homogeneous Neumann boundary conditions is solved and used to construct the superconducting-normal phase boundary for each sample. Vortex-antivortex (V-AV) configurations for each sample that accurately reflect the discrete symmetry of each sample boundary were found through the computational approach. These V-AV configurations are realized just within the phase boundary in the magnetic field-temperature phase diagram. Comparisons are made between the results obtained for the different sample shapes.