In this study, Bayesian optimization was combined with experiments to optimize the experimental conditions for the synthesis of LaFeO3 crystals with desired crystal properties using the flux growth method. The first set of candidates for the experimental conditions was determined based on the design of experiments. LaFeO3 crystals were synthesized under specific conditions, and their crystal properties were analyzed. Using the obtained data set, a Gaussian process regression model Y = f(X) was constructed with experimental conditions as explanatory variables X and crystal properties as objective variables Y. Candidate experimental conditions were then input into the model to predict crystal properties, and an acquisition function was calculated. Candidate experimental conditions with large values of the acquisition function were selected. Subsequently, LaFeO3 crystals were synthesized under these conditions, and their crystal properties analyzed. Construction of the Gaussian process regression model, selection of the next candidates for experimental conditions, and synthesis were repeated to synthesize LaFeO3 crystals with the desired crystal properties. Three case studies confirmed that LaFeO3 crystals with small crystallite sizes, calculated using X-ray diffraction (XRD) patterns and Scherrer's equation, can be synthesized within a small number of experiments. Furthermore, an index was prepared to evaluate the presence of impurities by excluding the LaFeO3-derived peaks from the XRD patterns of the products, which were used as Y in addition to crystallite size. A multiobjective optimization was conducted, and the suppression of impurity formation and reduction of crystallite size were achieved simultaneously.