The use of B-spline interpolation functions in the finite element method (FEM) is not a new subject. B-splines have been utilized in finite elements for many reasons. One reason is the higher continuity of derivatives and smoothness of B-splines. Another reason is the possibility of reducing the required number of degrees of freedom compared to a conventional finite element analysis. Furthermore, if B-splines are utilized to represent the geometry of a finite element model, interfacing a finite element analysis program with existing computer aided design programs (which make extensive use of B-splines) is possible. While B-splines have been used in finite element analysis due to the aforementioned goals, it is difficult to find resources that describe the process of implementing B-splines into an existing finite element framework. Therefore, it is necessary to document this methodology. This implementation should conform to the structure of conventional finite elements and only require exceptions in methodology where absolutely necessary. One goal is to implement B-spline interpolation functions in a finite element framework such that it appears very similar to conventional finite elements and is easily understandable by those with a finite element background. The use of B-spline functions in finite element analysis has been studied for advantages and disadvantages. Two-dimensional B-spline and standard FEM have been compared. This comparison has addressed the accuracy as well as the computational efficiency of B-spline FEM. Results show that for a given number of degrees of freedom, B-spline FEM can produce solutions with lower error than standard FEM. Furthermore, for a given solution time and total analysis time B-spline FEM will typically produce solutions with lower error than standard FEM. However, due to a more coupled system of equations and larger elemental stiffness matrix, B-spline FEM will take longer per degree of freedom for solution and assembly times than standard FEM. Three-dimensional B-spline FEM has also been validated by the comparison of a three-dimensional model with plane-strain boundary conditions to an equivalent two-dimensional model using plane strain conditions.