Although Bortfeld's analytical formula is useful for describing Bragg curves, measured data can deviate from the values predicted by the model. Thus, we sought to determine the parameters of a closed analytical expression of multiple Bragg curves for scanning proton pencil beams using a simultaneous optimization algorithm and to determine the minimum number of energies that need to be measured in treatment planning so that complete Bragg curves required by the treatment planning system (TPS) can be accurately predicted. We modified Bortfeld's original analytical expression of Bragg curves to accurately describe the dose deposition resulting from secondary particles. The parameters of the modified analytical expression were expressed as the parabolic cylinder function of the ranges of the proton pencil beams in water. Thirty-nine discrete Bragg curves were measured in our center using a PTW Bragg Peak chamber during acceptance and commission of the scanning beam proton delivery system. The coefficients of parabolic function were fitted by applying a simultaneous optimization algorithm to seven measured curves. The required Bragg curves for 45 energies in the TPS were calculated using our parameterized analytical expression. Finally, the 10 cm width of spread-out Bragg peaks (SOBPs) of beams with maximum energies of 221.8 and 121.2 MeV were then calculated in the TPS and compared with measured data. Compared with Bortfeld's original formula, our modified formula improved fitting of the measured depth dose curves at depths around three-quarters of the maximum range and in the beam entrance region. The parabolic function described the relationship between the parameters of the analytic expression of different energies. The predicted Bragg curves based on the parameters fitted using the seven measured curves accurately described the Bragg curves of proton pencil beams of 45 energies configured in our TPS. When we used the calculated Bragg curves as the input to TPS, the standard deviations of the measured and calculated data points along the 10 cm SOBPs created with proton pencil beams with maximum energies of 221.8 and 121.2 MeV were 1.19% and 1.18%, respectively, using curves predicted by the algorithm generated from the seven measured curves. Our method would be a valuable tool to analyze measured Bragg curves without the need for time-consuming measurements and correctly describe multiple Bragg curves using a closed analytical expression.