We consider the electroweak production of multiple Higgs bosons at the LHC in the two-Higgs doublet model (2HDM). We show that in realistic models the electroweak cross sections for production of 2 or 3 Higgs bosons can dominate over QCD production. This is because in a gg-initiated process, production via a resonant Z in the s-channel is prohibited according to the Landau-Yang theorem, which is not the case for a qq-initiated process. We discuss the issue of accessing the 3-particle couplings with combinations of Higgs bosons and vector bosons, as well as Higgs-fermion couplings, using electroweak production modes in final states with multiple Higgs bosons and/or vector bosons. We therefore explore regions of parameter space of the Type-I and Type-II models that are consistent with constraints from collider searches, b-physics and EW precision data and the data on the Higgs boson at LHC. We compute cross sections for Higgs pair production qq --> h_i h_j and associated production with a gauge boson, qq --> h_i V_j, where h_i is any of the 2HDM scalars. We also consider production of 3 and 4 Higgs bosons, and we compute the corresponding gg initiated cross sections. In some interesting cases the qq initiated cross sections can be as large as hundreds of femtobarns. We also consider the special case when the Higgs boson discovered at the LHC is the heavier CP-even Higgs boson of the 2HDM, and when the sum of the masses of the lighter CP-even scalar and the CP-odd scalar is smaller than the Z boson mass. We find that there is a region of parameter space in the Type-I model where this is possible, and study the phenomenology of such parameter points.