This paper is an extended version of the ICCS conference paper (Kapturczak et al., 2020, [1]) and presents a way to improve the boundary shape modeling process in solving boundary value problems in elasticity. The inclusion of NURBS curves into the mathematical formalism of the parametric integral equations system method (PIES) is proposed. The advantages of such an application are widely discussed. Recently, the Bezier curves, mainly the cubic curves (of third-degree), were used. The segments of the boundary shape were modeled by such curves (with ensuring continuity at the connection points). Using NURBS curves, the boundary shape can be modeled with only one curve. So, continuity is automatically ensured. Additionally, the second degree NURBS curve is enough to obtain the shape with high accuracy (better than cubic Bezier curves). The NURBS curve is defined by points, their weights, and the knots vector. Such parameters significantly improve the shape modification process, which can directly improve e.g. the shape identification process. The examples of shape modification using such parameters are presented. The boundary shapes of the examples (even defined by both linear and curvilinear segments) can be defined using only one closed NURBS curve. The impact of modeling accuracy on the final PIES solutions is examined on examples described by the Navier–Lamé equations. To improve calculations, the PIES method using NURBS curves was implemented as a computer program. Then, it was decided to verify the accuracy of the obtained solutions. For comparison, the solutions were also obtained using analytical solutions, boundary element method, and PIES method (with the Bezier curves). An improvement in the boundary shape modeling was noticed. It significantly affects the accuracy of solutions. As a result, the consumption of computer resources was reduced, while the process of boundary shape modeling and the accuracy of the obtained results were improved.