A general theory of piecewise multiharmonic splines is constructed for a class of fractals (post-critically finite) that includes the familiar Sierpinski gasket, based on Kigami's theory of Laplacians on these fractals. The spline spaces are the analogues of the spaces of piecewise Cj polynomials of degree 2j + 1 on an interval, with nodes at dyadic rational points. We give explicit algorithms for effectively computing multiharmonic functions (solutions of Δj+1u = 0) and for constructing bases for the spline spaces (for general fractals we need to assume that j is odd), and also for computing inner products of these functions. This enables us to give a finite element method for the approximate solution of fractal differential equations. We give the analogue of Simpson's method for numerical integration on the Sierpinski gasket. We use splines to approximate functions vanishing on the boundary by functions vanishing in a neighbourhood of the boundary.