What is the difference between cubic interpolation and cubic "Spline" interpolation?. How to use it for upsampling purpose?
The difference between cubic interpolation as described in your question and cubic spline interpolation is that in cubic interpolation you use 4 data points to compute the polynomial. There are no constraints on the derivatives. Cubic spline interpolation computes a third order polynomial only from two data points with the additional constraint that the first and second derivative at the interpolation points are continuous. So if you have 4 points, then you compute 3 different polynomials (between points 1-2, 2-3, and 3-4), and these polynomials are smoothly connected in the sense that their first and second derivatives are equal at the given data points. This is not the case with standard polynomial interpolation.

But if you want an efficient implementation I wouldn't necessarily use polynomial or spline interpolation but I would go for a standard polyphase implementation of a low-pass interpolation filter. If $L$ is the upsampling factor you get $L$ parallel polyphase filters operating at the lower sampling rate, followed by $L$ upsamplers, the results of which are added. The low-pass filter can be designed using the window method (e.g. a Kaiser window). There is a lot of material available on the polyphase implementation of interpolators. This Master's thesis is very accessible. Just have a look at single-stage polyphase implementations using linear-phase FIR filters and ignore the rest.

EDIT:

If you want or need to use cubic interpolation, here's what you need to do. First of all, some coefficients in your equations are wrong, that's why also the result is not correct. These would be the correct equations:


    V(2) = a + 2b + 4c + 8d
    V(3) = a + 3b + 9c + 27d
    V(4) = a + 4b + 16c + 64d
    V(5) = a + 5b + 25c + 125d


However, you don't want to do it this way because if you always use the sample index as the $x$-coordinate you'll end up with a very large range for the coefficients of the system of linear equations (because you need to take the third power), and your coefficient matrix will sooner or later become ill-conditioned (asking for numerical problems). Second, and even more importantly, you would get a different system of equations for every new segment of the input signal, which would be very inefficient. You can do the interpolation by solving a $4\times 4$ system of linear equations with the same coefficient matrix each time. Just the left-hand side vector (containing the values of the given data points, V in your example) will be different for each segment.

You compute a cubic polynomial for each segment (i.e. for each range between two given data points) by considering the data points defining the segment and the two adjacent data points, just as in your example. If you define the cubic polynomial as

$$P(x)=a_0+a_1x+a_2x^2+a_3x^3$$

and if you define the $x$-axis values of the 4 current input data points as $-1$, $0$, $1$, and $2$, you get the following system of equations:

$$\begin{align}y_0&amp;=a_0-a_1+a_2-a_3\\
y_1&amp;=a_0\\
y_2&amp;=a_0+a_1+a_2+a_3\\
y_3&amp;=a_0+2a_1+4a_2+8a_3\end{align}$$

where $y_i$, $k=0,1,\ldots, 3$ are the current data points, and the current segment is between the points $y_1$ and $y_2$. The system above is readily solved as

$$\begin{align}
a_0&amp;=y_1\\
a_1&amp;=-\frac13 y_0 - \frac12 y_1 + y_2 - \frac16 y_3\\
a_2&amp;=\frac12  (y_0 + y_2) - y_1\\
a_3&amp;=\frac12  \left[(y_1 - y_2) + \frac13 (y_3 - y_0) \right]\\
\end{align}$$

So you just need to evaluate the above 4 expressions for the polynomial coefficients for each set of 4 data values. As soon as you have the coefficients, you can compute the interpolated output values in the current segment by evaluating the polynomial at $x_i=i/L$, $i=1,2,\ldots,L-1$, where $L$ is the upsampling factor ($L=10$ in your case). Then you forget the oldest input value, read the new input value, and you have a new set of 4 values from which you compute a new polynomial for the current segment, etc.

Here is a little Matlab/Octave script implementing this method:


x = randn(20,1);        % some input data
L = 10;                 % upsampling factor

% append points at both ends for first and last segment
N = length(x);
x2 = [x(1);x(:);x(N)];
Nseg = N-1;         % # segments

% matrix to compute polynomial coefficients from the current 4 data points: a = A*y
A = [0 1 0 0;-1/3 -.5 1 -1/6;.5 -1 .5 0;-1/6 .5 -.5 1/6];

% x-axis values of interpolated values in the current segment
mu = (0:L-1)/L;

% allocate output vector
Ny = Nseg*L + 1;
y = zeros(Ny,1);

% loop over all segments
for i = 1:Nseg,
    xt = x2(i:i+3);
    a = A*xt;
    yt = a(1) + a(2)*mu + a(3)*mu.^2 + a(4)*mu.^3;
    Iy = (i-1)*L + (1:L);
    y(Iy) = yt;
end

y(Ny) = x(N);

% plot result
tx = 0:N-1;
ty = (0:Ny-1)*(N-1)/(Ny-1);
plot(ty,y,'r-'), hold on
stem(ty,y,'r'), hold on
stem(tx,x), hold off