Owing to the significant increase in the volume of world trade, mega-container vessels are being used to meet transportation demands. As the size of the vessels increases, the loading sequence of containers onto the vessels presents an important challenge for planners. In this study, we consider the container stowage planning problem with stability constraints (e.g. shear force, bending moment, trim) and develop a mixed integer linear programming (MILP) formulation which generates load plans by minimizing total cost associated with the over-stows and trimming moments. Our study adopts a holistic perspective which encompasses several real-world features such as different container specifications, a round-robin tour of multiple ports, technical limitations related to stack weight, stress, and ballast tanks. We also propose a two-stage heuristic solution methodology that employs an integer programming (IP) formulation and then a swapping heuristic (SH) algorithm. This approach first acquires a lower bound on the total over-stow cost with the IP model, thereby creating an initial bay plan. Then, it applies the SH algorithm to this initial bay plan to minimize cost resulting from trimming moments. The efficiency of the MILP formulation and heuristic algorithm is investigated through numerical examples. The results have shown that the heuristic has greatly improved the solution times as well as the size of the solvable problems compared to the MILP formulation. In particular, the two-stage heuristic can solve all size problem instances within an average optimality gap of 0-25% in less than 8 minutes, whereas the MILP can only achieve an approximate optimality gap of 55-80% in 2 hours.