In this thesis we implement and employ a simplex spline basis developed by Cohen, Lyche, and Riesenfeld for the space of C1 quadratic splines on the Powell–Sabin 12-split of a triangle. A matrix recursion is used for evaluation and differentiation, where the need for standard Bernstein– Bézier techniques are avoided. A brief account of the construction of multivariate splines and the theoretical foundation underlying the finite element method is given. Subsequently, using an explicit conversion to the Hermite nodal basis known from the finite element method, we solve the biharmonic equation using conforming methods. Interpolation estimates are derived, and numerical results are seen to comply.