This thesis presents the search for heavy neutral Higgs bosons, produced in association with b-quarks and decaying into a b-quark pair in proton-proton collisions at a centre-of-mass energy of 13 TeV at the LHC with the CMS detector. The process is chosen, given that a considerable enhancement of the coupling to b-quarks is possible in theories beyond the Standard Model, such as the Minimal Supersymmetric extension of the Standard Model (MSSM) and the Two Higgs Doublet Model (2HDM). A new boson would manifest itself as an excess in the invariant mass distribution of the two leading b-tagged jets. The final state is rather challenging, requiring dedicated triggers to control the large rate of multijet backgrounds. The data used in this thesis was collected during 2017 and 2018 and combined with a previous search in the same channel with 2016 data to provide results using the full Run 2 dataset, accounting for up to 127~/ifb of data. The search is performed for several mass hypotheses in two categories. In one, it is conducted in a fully-hadronic category, which is sensitive only to masses above 300 GeV and limited by trigger rates. Sensitivity to lower masses is recovered in the semi-leptonic category, in which a muon is required within a final state jet, reducing the trigger rates with low transverse-momenta jets. Special attention is given to the background modelling, for which a data-driven approach has been used, and efforts have been made to reduce the uncertainties. In addition, a new modified event selection yields significant enhancements in the signal efficiency and overall significance compared to previous analyses. These improvements enabled the exploration of masses up to 1.8 TeV, a significant extension compared to existing results. The sensitivity towards lower masses was extended to 125 GeV with the SL category. No significant deviation with respect to the background-only hypothesis was observed, and upper limits were set on the cross-section times branching ratio of the scrutinized process at 95/% confidence level. These limits are translated into exclusion limits in the parameter space of MSSM and 2HDM benchmark scenarios. These are the most stringent limits to date on this process.