We present a search for neutral Higgs bosons ϕ decaying into bb, produced in association with b quarks in pp collisions. This process could be observable in supersymmetric models with high values of tanβ. The event sample corresponds to 2.6  fb−1 of integrated luminosity collected with the CDF II detector at the Fermilab Tevatron collider. We search for an enhancement in the mass of the two leading jets in events with three jets identified as coming from b quarks using a displaced vertex algorithm. A data-driven procedure is used to estimate the dijet mass spectrum of the nonresonant multijet background. The contributions of backgrounds and a possible Higgs boson signal are determined by a two-dimensional fit of the data, using the dijet mass together with an additional variable which is sensitive to the flavor composition of the three tagged jets. We set mass-dependent limits on σ(pp→ϕb)×B(ϕ→bb) which are applicable for a narrow scalar particle ϕ produced in association with b quarks. We also set limits on tanβ in supersymmetric Higgs models including the effects of the Higgs boson width.