The paper studies the problem of construction of optimal interpolation formulas with derivative in the Sobolev space L (m) 2 (0, 1). Here the interpolation formula consists of the linear combination of values of the function at nodes and values of the first derivative of that function at the end points of the interval [0, 1]. For any function of the space L (m) 2 (0, 1) the error of the interpolation formulas is estimated by the norm of the error functional in the conjugate space L (m) * 2 (0, 1). For this, the norm of the error functional is calculated. Further, in order to find the minimum of the norm of the error functional, the Lagrange method is applied and the system of linear equations for coefficients of optimal interpolation formulas is obtained. It is shown that the order of convergence of the obtained optimal interpolation formulas in the space L (m) 2 (0, 1) is O(h m). In order to solve the obtained system it is suggested to use the Sobolev method which is based on the discrete analog of the differential operator d 2m / dx 2m. Using this method in the cases m = 2 and m = 3 the optimal interpolation formulas are constructed. It is proved that the order of convergence of the optimal interpolation formula in the case m = 2 for functions of the space C 4 (0, 1) is O(h 4) while for functions of the space L (2) 2 (0, 1) is O(h 2). Finally, some numerical results are presented,