  Cubic spline interpolation on Euclidean space is a standard topic in
numerical analysis, with countless applications in science and technology. In
several emerging fields, for example computer vision and quantum control, there
is a growing need for spline interpolation on curved, non-Euclidean space. The
generalization of cubic splines to manifolds is not self-evident, with several
distinct approaches. One possibility is to mimic the acceleration minimizing
property, which leads to Riemannian cubics. This, however, require the solution
of a coupled set of non-linear boundary value problems that cannot be
integrated explicitly, even if formulae for geodesics are available. Another
possibility is to mimic De~Casteljau's algorithm, which leads to generalized
B\'ezier curves. To construct C2-splines from such curves is a complicated
non-linear problem, until now lacking numerical methods. Here we provide an
iterative algorithm for C2-splines on Riemannian symmetric spaces, and we prove
convergence of linear order. In terms of numerical tractability and
computational efficiency, the new method surpasses those based on Riemannian
cubics. Each iteration is parallel, thus suitable for multi-core
implementation. We demonstrate the algorithm for three geometries of interest:
the $n$-sphere, complex projective space, and the real Grassmannian.
