In this paper, Chebyshev series and rigorous numerics are combined to compute solutions of the Euler-Lagrange equations for the one-dimensional Ginzburg-Landau model of superconductivity. The idea is to recast solutions as fixed points of a Newton-like operator defined on a Banach space of rapidly decaying Chebyshev coefficients. Analytic estimates, the radii polynomials and the contraction mapping theorem are combined to show existence of solutions near numerical approximations. Coexistence of as many as seven nontrivial solutions is proved.