We investigated the accuracy and range of validity of the perturbative model for the two-point (2PCF) and three-point (3PCF) correlation functions in real space in view of the forthcoming analysis of the Euclid mission spectroscopic sample. We took advantage of clustering measurements from four snapshots of the Flagship I N-body simulations at z = {0.9,1.2,1.5,1.8}, which mimic the expected galaxy population in the ideal case, i.e. in the absence of observational effects such as purity and completeness. For the 3PCF we considered all available triangular configurations given a minimal separation (rmin). We first assessed the model performance by fixing the cosmological parameters and evaluating the goodness of fit provided by the perturbative bias expansion in the joint analysis of the two statistics, finding an overall agreement with the data down to separations of 20 h−1 Mpc. Subsequently, we built on the state-of-the-art analysis and extended it to include the dependence on three cosmological parameters: the amplitude of scalar perturbations (As), the matter density (ωcdm), and the Hubble parameter (h). To achieve this goal, we developed an emulator capable of generating fast and robust modelling predictions for the two summary statistics, which thus enables an efficient sampling of the joint likelihood function. We therefore present the first joint full-shape analysis of the real-space 2PCF and 3PCF, testing the consistency and constraining power of the perturbative model across both probes and assessing its performance in a combined likelihood framework. We explored possible systematic uncertainties induced by the perturbative model at small scales, finding an optimal scale cut of rmin = 30 h−1 Mpc for the 3PCF when imposing an additional limitation on the nearly isosceles triangular configurations included in the data vector. This work is part of a series of papers in which we validate theoretical models for galaxy clustering measurements in preparation for the Euclid mission.