Nowadays, the finite element method (FEM) is the most known and the most commonly used method for solving boundary value problems defined exactly (without uncertainty). This method is based on the division of the domain into finite elements. These elements are connected to each other and produces a mesh. The second method, more recent, but also often used, is the boundary element method (BEM). This method is based on the division of the boundary only on the so-called boundary elements. This solution significantly reduces the amount of data required to solve the problem. Therefore, the amount of computer resources, needed to obtain a solution, is also significantly reduced. However, both FEM and BEM require a large number of points to define the shape of the boundary. It was the main reason to look for new methods. This is how the parametric integral equations system (PIES) appeared. In this paper, the application of interval arithmetic for modeling uncertainly defined boundary conditions (in solving boundary value problems) is presented. In the literature, for solving such problems, different modifications of mentioned FEM and BEM methods appeared. Many different ways of modeling uncertainty have been used. However, the most important problem is to reduce the number of calculations. This is the main reason to choose interval arithmetic. Direct application of classical or directed interval arithmetic, into interval modifications of FEM and BEM methods, leads to significant overestimation. Therefore, it was decided to propose interval modification of the parametric integral equations system (PIES) method. Recently, the PIES method has been thoroughly tested on examples of boundary problems, where input data were defined without uncertainty. During tests, the number of input data (necessary to define the problem) was always smaller (comparing with well-known methods). This also allows reducing the number of equations in the obtained system of equations. This gives a significant advantage (over other methods) in modeling uncertainty of input data. The application of the interval PIES method resulted in a significant reduction of overestimation. However, obtained solutions (by direct application of interval arithmetic known from literature) are not exactly in line with expectations. Therefore, the modification of directed interval arithmetic has been proposed. It has been applied for calculations in the interval PIES method. The proposed method has been tested on examples of boundary value problems (modeled by Laplace's equation) with uncertainly defined boundary conditions. The uncertainty of boundary conditions has been defined as interval constant boundary conditions as well as interval linear function. Obtained solutions have been compared with the solutions of exactly defined (without uncertainty) PIES.