In this work, a numerical–experimental study of the interlaminar zone for an unidirectional glass fiber reinforced epoxy composite is carried out in order to predict the load–displacement curves of a double cantilever beam test. First, an experimental mechanical characterization of the laminated composite was made through quasi-static in-plane tensile and bending tests and out-of-plane delamination tests (i.e. double cantilever beam tests or opening mode I). The main results have shown that the elastic module in the fiber direction is E1 = 32.1 GPa and the fracture process is characterized by a critical energy release rate G IC = 1466 N/m. In order to be able to predict delamination using finite element analysis, a bilinear cohesive zone model is adopted. This law has two parameters—the initial stiffness [Formula: see text] and the critical traction [Formula: see text]—to be fit for reducing the distance between numerical and experimental double cantilever beam load–displacement curves. That means an optimization problem, which is here solved by proposing a very simple and cheap three-step procedure, avoiding expensive three-dimensional simulations: (i) the numerical double cantilever beam test is coded in the OCTAVE software using one-dimensional beam elements and the bilinear cohesive zone model; (ii) several one-dimensional simulations are performed varying at the cohesive law's parameters in order to build the objective function; and (iii) a genetic algorithm from the Scilab optimization toolbox is then used to determine the interface's parameters minimizing the objective function. This method has proposed [Formula: see text] = 2.74 × 1014 N/m3 and [Formula: see text] = 1.99 MPa as the best parameters fitting the curves. Finally, the identified bilinear cohesive law's parameters are used to carry out a full three-dimensional finite element analysis simulation in ANSYS, showing the same response than the simplified beam code, but with the possibility to obtain more accurate information about the fracture process (e.g. the shape of the process zone) or to use more complicated geometries or boundary conditions. In this work, we show the ability to predict the material behavior through different methods, using the best characteristics of each one.