C++ code implements http://en.wikipedia.org/w/index.php?title=Spline_%28mathematics%29&oldid=288288033#Algorithm_for_computing_natural_cubic_splines

Plot-Splint.pyw

import numpy as np
from pylab import *

raw = np.array([
	[0, 0.0],
	[0.15708, 0.156313],
	[0.314159, 0.308879],
	[0.471239, 0.453952],
	[0.628319, 0.587785],
	[0.785398, 0.706871],
	[0.942478, 0.808655],
	[1.09956, 0.890822],
	[1.25664, 0.951057],
	[1.41372, 0.987429],
	[1.5708, 0.999553],
	[1.72788, 0.987429],
	[1.88496, 0.951057],
	[2.04204, 0.890822],
	[2.19911, 0.808655],
	[2.35619, 0.706871],
	[2.51327, 0.587785],
	[2.67035, 0.453952],
	[2.82743, 0.308879],
	[2.98451, 0.156313],
	[3.14159, 5.35898e-008],
	[3.29867, -0.156313],
	[3.45575, -0.308879],
	[3.61283, -0.453952],
	[3.76991, -0.587785],
	[3.92699, -0.706871],
	[4.08407, -0.808655],
	[4.24115, -0.890822],
	[4.39823, -0.951056],
	[4.55531, -0.987429],
	[4.71239, -0.999553],
	[4.86947, -0.987429],
	[5.02655, -0.951057],
	[5.18363, -0.890822],
	[5.34071, -0.808655],
	[5.49779, -0.706871],
	[5.65487, -0.587785],
	[5.81195, -0.453952],
	[5.96903, -0.308879],
	[6.12611, -0.156313],
	[6.28319, -1.0718e-007],
	]);

import numpy as np
import matplotlib.pyplot as plt
#raw = np.array([[0,0],[1,0.841471],[2,0.909297],[3,0.14112],[4,-0.756802],[5,-0.958924],[6,-0.279415],[7,0.656987],[8,0.989358],[9,0.412118]]);
x = raw[:,0]
y = raw[:,1]
f_x = np.linspace(x[0], x[-1], 100)
f_sin = np.sin(f_x)

ax = gca()
ax.spines['right'].set_color('none') # discard right spines
ax.spines['top'].set_color('none') # discard top spines

# move top spines
ax.xaxis.set_ticks_position('bottom')
ax.spines['bottom'].set_position(('data',0))
# move left spines
ax.yaxis.set_ticks_position('left')
ax.spines['left'].set_position(('data',0))
axis([0, 2 * pi + 0.5, -1.1, 1.1])
xticks([-pi/4, 0, pi/2, pi, pi * 2],
    [r'$-\pi/4$', r'$0$', r'$+\pi/2$', r'$+\pi$', r'$2\pi$'])

yticks([-1, 0, +1], [r'$-1$', r'$0$', r'$+1$'])
plt.plot(x, y, "o", f_x, f_sin, "--", x, y, "-")
plt.legend(['data', "sin(x)", 'spline'], loc='best')
plt.show()

spline.cpp

     // http://en.wikipedia.org/w/index.php?title=Spline_%28mathematics%29&oldid=288288033#Algorithm_for_computing_natural_cubic_splines
    #include <iostream>
    #include <vector>
    #include <algorithm>
    #include <cmath>
    #include <cstdio>
    using namespace std;

    typedef vector<double> vec;

    struct SplineSet {
        double a;
        double b;
        double c;
        double d;
        double x;
    };

    // Solve Spline Si
    vector<SplineSet> spline(vec &x, vec &y)
    {
        int n = x.size() - 1;
        vec a;
        a.insert(a.begin(), y.begin(), y.end());
        vec b(n);
        vec d(n);
        vec h;

        for(int i = 0; i < n; ++i)
            h.push_back(x[i+1]-x[i]);

        vec alpha;
        for(int i = 0; i < n; ++i)
            alpha.push_back( 3*(a[i+1]-a[i])/h[i] - 3*(a[i]-a[i-1])/h[i-1]  );

        vec c(n+1);
        vec l(n+1);
        vec mu(n+1);
        vec z(n+1);
        l[0] = 1;
        mu[0] = 0;
        z[0] = 0;

        for(int i = 1; i < n; ++i)
        {
            l[i] = 2 *(x[i+1]-x[i-1])-h[i-1]*mu[i-1];
            mu[i] = h[i]/l[i];
            z[i] = (alpha[i]-h[i-1]*z[i-1])/l[i];
        }

        l[n] = 1;
        z[n] = 0;
        c[n] = 0;

        for(int j = n-1; j >= 0; --j)
        {
            c[j] = z [j] - mu[j] * c[j+1];
            b[j] = (a[j+1]-a[j])/h[j]-h[j]*(c[j+1]+2*c[j])/3;
            d[j] = (c[j+1]-c[j])/3/h[j];
        }

        vector<SplineSet> output_set(n);
        for(int i = 0; i < n; ++i)
        {
            output_set[i].a = a[i];
            output_set[i].b = b[i];
            output_set[i].c = c[i];
            output_set[i].d = d[i];
            output_set[i].x = x[i];
        }
        return output_set;
    }


    // calculate the Si(x) = a + b(x - xi) + c(x - xi)^2 + d(x- xi)^3
    // optimized to 
    double interpolate(SplineSet const &cs, double x, double xi)
    {
        double delta = x - xi;
        return cs.a + (cs.b + (cs.c + cs.d * delta) * delta) * delta;
    }

int main()
{
    vec x(11);
    vec y(11);
    const double pi = 3.1415926;
    //  0 ~ 2pi;
    for(int i = 0; i < x.size(); ++i) {
    	x[i] = 2 * pi * i / 10.0;
 	    y[i] = sin(x[i]);
 	    cout << x[i] << ", " << y[i] << endl;
    }
  	cout << endl;
    vector<SplineSet> cs = spline(x, y);
    for(int i = 0; i < cs.size(); ++i)
        cout << cs[i].d << "\t" << cs[i].c << "\t" << cs[i].b << "\t" << cs[i].a << endl;

    cout << endl;
    for(int i = 1; i < x.size(); ++i) {
        double cx = x[i - 1];
        double cy = y[i - 1];
        double xt = cx;
        double yt = cy;
        for (int j = 0; j < 4; ++j)
        {
           xt += pi / 20.0; 
           yt += pi / 20.0;
    	   cout << xt << ", "
                << interpolate(cs[i - 1], yt, cy) << endl;
        }
    }
    return getchar();
}