Taylor series of a polynomial Given a polynomial $y=C_0+C_1 x+C_2 x^2+C_3 x^3 + \ldots$ of some order $N$, I can easily calculate the polynomial of reduced order $M$ by taking only the first $M+1$ terms. This is equivalent to doing a Taylor series expansion with $M<=N$ around $x=0$.
But what if I want to take the Taylor series expansion around a different point $x_c$. In the end, I want the polynomial coefficients of $y_2=K_0+K_1 x + K_2 x^2 + K_3 x^3 + \ldots$ which represents the Taylor's expansion of $y$ around point $x_c$ such that $y(x_c)=y_2(x_c)$ including the first $M$ derivatives.
So given the coefficients $C_i$ with $i=0 \ldots N$, and a location $x_c$ I want to calculate the coefficients $K_j$ with $j=0 \ldots M$.
Example
Given $y=C_0+C_1 x+C_2 x^2$ ( $N=2$ ) then the tangent line ($M=1$) through $x_c$ is
$$ y_2 = (C_0-C_2 x_c^2) + (C_1+2 C_2 x_c) x $$
or $K_0 = C_0-C_2 x_c^2$, and $K_1 =C_1+2 C_2 x_c$
There must be a way to construct a ($M+1$ by $N+1$ ) matrix that transforms the coefficients $C_i$ into $K_j$. For the above example this matrix is
$$ \begin{bmatrix}K_{0}\\
K_{1}\end{bmatrix}=\begin{bmatrix}1 & 0 & -x_{c}^{2}\\
0 & 1 & 2\, x_{c}\end{bmatrix}\begin{bmatrix}C_{0}\\
C_{1}\\
C_{2}\end{bmatrix} $$
Example #2
The reduction of a $5$-th order polynomial to a $3$-rd order around $x_c$ is
$$ \begin{bmatrix}K_{0}\\
K_{1}\\
K_{2}\\
K_{3}\end{bmatrix}=\left[\begin{array}{cccc|cc}
1 &  &  &  & -x_{c}^{4} & -4\, x_{c}^{5}\\
 & 1 &  &  & 4\, x_{c}^{3} & 15\, x_{c}^{4}\\
 &  & 1 &  & -6\, x_{c}^{2} & -20\, x_{c}^{3}\\
 &  &  & 1 & 4\, x_{c} & 10\, x_{c}^{2}\end{array}\right]\begin{bmatrix}C_{0}\\
C_{1}\\
C_{2}\\
C_{3}\\
C_{4}\\
C_{5}\end{bmatrix} $$
which is a block matrix, and not an upper diagonal one as some of the answers have indicated.
 A: Here is a Mathematica routine that is (more or less) an efficient way of performing Arturo's proposal (I assume the array of coefficients cofs is arranged with constant term first, i.e. $p(x)=\sum\limits_{k=0}^n$cofs[[k + 1]]$x^k$):
polxpd[cofs_?VectorQ, h_, d_] := Module[{n = Length[cofs] - 1, df},
    df = PadRight[{Last[cofs]}, d + 1];
    Do[
       Do[
          df[[j]] = df[[j - 1]] + h df[[j]],
          {j, Min[d, n - k + 1] + 1, 2, -1}];
       df[[1]] = cofs[[k]] + h df[[1]],
       {k, n, 1, -1}];
    Do[
       Do[
          df[[k]] -= h df[[k + 1]],
          {k, d, j, -1}],
       {j, d}];
    df]

Let's try it out:
polxpd[{c[0], c[1], c[2]}, h, 1] // FullSimplify
{c[0] - h^2*c[2], c[1] + 2*h*c[2]}

polxpd[c /@ Range[0, 5], h, 3] // FullSimplify
{c[0] - h^4*(c[4] + 4*h*c[5]), c[1] + h^3*(4*c[4] + 15*h*c[5]), 
 c[2] - 2*h^2*(3*c[4] + 10*h*c[5]), c[3] + 2*h*(2*c[4] + 5*h*c[5])}

Now, Arturo gave the linear-algebraic interpretation of this conversion; I'll look at this from the algorithmic point of view:
For instance, see this (modified) snippet:
n = Length[cofs] - 1;
df = {Last[cofs]};
Do[
   df[[1]] = cofs[[k]] + x df[[1]],
   {k, n, 1, -1}];

This is nothing more than the Horner scheme (alias "synthetic division") for evaluating the polynomial at x. What is not so well known is that the Horner scheme can be hijacked so that it computes derivatives as well as polynomial values. We can "differentiate" the previous code snippet like so (i.e., automatic differentation):
n = Length[cofs] - 1;
df = {Last[cofs], 0};
Do[
   df[[2]] = df[[1]] + x df[[2]];
   df[[1]] = cofs[[k]] + x df[[1]],
  {k, n, 1, -1}];

where the rule is $\frac{\mathrm d}{\mathrm dx}$df[[j]]$=$df[[j+1]]. "Differentiating" the line df = {Last[cofs]} (the leading coefficient of the polynomial) requires appending a 0 (the derivative of a constant is $0$); "differentiating" the evaluation line df[[1]] = cofs[[k]] + x df[[1]] gives df[[2]] = df[[1]] + x df[[2]] (use the product rule, and the fact that $\frac{\mathrm d}{\mathrm dx}$cofs[[k]]$=0$). Continuing inductively (and replacing the x with h), we obtain the first double loop of polxpd[].
Actually, the contents of df after the first double loop are the "scaled derivatives"; that is, df[[1]]$=p(h)$, df[[2]]$=p^\prime(h)$, df[[3]]$=\frac{p^{\prime\prime}(h)}{2!}$, ... and so on.
What the second double loop accomplishes is the "shifting" of the polynomial by -h; this is in fact synthetic division applied repeatedly to the coefficients output by the first double loop, as mentioned here.
