Higher-order integral equation methods have been investigated. The study has focused on improving the accuracy and efficiency of the Method of Moments (MoM) applied to electromagnetic problems. A new set of hierarchical Legendre basis functions of arbitrary order is developed. The new basis functions provide much better accuracy than low-order basis functions for a fixed number of functions. The hierarchical feature of the new basis functions combines the advantages of low- and higher-order basis functions into a single flexible set of functions. In comparison to existing higher-order basis functions, the new basis functions result in a much lower matrix condition number which is accomplished by focusing on the orthogonality of the basis functions. Numerical results confirm that higher-order convergence and very favorable condition numbers are achieved. The low matrix condition number obtained with the new basis functions enables an efficient iterative solution of higher-order MoM systems. Iterative solution methods incorporate a matrix preconditioner and four preconditioners are presented here; two of these are found in existing works and the other two are adaptations of existing preconditioners to the higher-order case. Numerical experiments verify that the combination of the new basis functions, a good preconditioner, and a robust iterative algorithm lead to fast convergence even for very high-order basis functions, e.g. 10th order. The new Legendre basis functions are further developed to incorporate edge singularities. Previous works in this area have shown the necessity of singular basis functions and three existing formulations are adapted to the new Legendre basis functions. Numerical experiments show that one formulation is the most accurate as well as the fastest to compute. Using this formulation, the accuracy of surface currents in the vicinity of edges is greatly improved and smaller improvements are obtained for the far field. The hybrid Physical Optics (PO) - MoM is an efficient technique for treating objects that are too large in terms of wavelengths to be analyzed with MoM. The existing hybrid technique employs low-order basis functions and flat patches. This technique is extended here to the case of higher-order hierarchical Legendre basis functions and curved patches. The required memory and the computation time of the new higherorder hybrid PO-MoM are typically reduced by a factor of 10 in comparison to the existing technique. The hybrid technique includes the coupling between the MoM and PO regions and numerical results are presented to illustrate the accuracy. The hierarchical feature of the new higher-order Legendre basis functions allows a flexible selection of the expansion order. This is demonstrated by numerical examples involving patches of non-uniform size. The expansion order is adapted to the electrical size of each patch which results in a very efficient solution.