We describe explicitly each stage of a numerically stable algorithm for calculating with exponential tension B-splines with non-uniform choice of tension parameters. These splines are piecewisely in the kernel of D 2 (D 2 − p 2 ), where D stands for ordinary derivative, defined on arbitrary meshes, with a different choice of the tension parameter p on each interval. The algorithm provides values of the associated B-splines and their generalized and ordinary derivatives by performing positive linear combinations of positive quantities, described as lower-order exponential tension splines. We show that nothing else but the knot insertion algorithm and good approximation of a few elementary functions is needed to achieve machine accuracy. The underlying theory is that of splines based on Chebyshev canonical systems which are not smooth enough to be ECC-systems. First, by de Boor algorithm we construct exponential tension spline of class C 1 , and then we use quasiOslo type algorithms to evaluate classical non-uniform C 2 tension exponential splines.