Plotting Partial Sums of Fourier Series

In the definition of s you're summing from k==0. Since the summand has a term 1/k this gives a divide-by-zero error when calculating the partial sums. The sum should in fact start from k==1 (the zeroth coefficient is taken care of by the constant term in front of the sum). The first few approximations then look like

s[n_, x_] := 
 8/4 + 3/(9 π) Sum[(6 (-1)^k)/(k π) Cos[(k π x)/
        2] + (16 (-1)^k + 13)/(π k) Sin[(k π x)/2], {k, 1, n}]

partialsums = Table[s[n, x], {n, 1, 5}];
Plot[Evaluate[partialsums], {x, -4, 4}]

Mathematica graphics

To compare this with f we plot the s[5,x] and f in the same plot:

Plot[{s[5, x], f[x]}, {x, -2, 2}]

Mathematica graphics

which doesn't seem right to me, so I suspect you made a mistake somewhere in calculating the coefficients.

Calculating coefficients by hand

We could use FourierSeries to calculate the partial sums, but this is very slow, and doesn't give produce the general equation for the coefficients. Therefore it's better to calculate the coefficients by hand. The $n$-th coefficient can be calculated according to

coeff[0] = 1/4 Integrate[f[x], {x, -2, 2}];
coeff[n_] = 1/4 Integrate[f[x] Exp[I Pi n x/2], {x, -2, 2}]
(1/(2 n^4 π^4))E^(-I n π) (-48 + 48 E^(I n π) - 
 48 I n π + 28 n^2 π^2 - 6 E^(I n π) n^2 π^2 + 
 2 E^(2 I n π) n^2 π^2 + 12 I n^3 π^3 - 
 I E^(I n π) n^3 π^3 - I E^(2 I n π) n^3 π^3)

Then the partial sums are given by

series[m_, x_] := Sum[Exp[-I Pi n x/2] coeff[n], {n, -m, m}]

Plotting the first few approximations:

Plot[Evaluate[Table[series[j, x], {j, 0, 5}]], {x, -6, 6}]

Mathematica graphics

To see how this compares with the original function f:

Plot[Evaluate[{series[5, x], f[Mod[x, 4, -2]]}], {x, -4, 4}]

Mathematica graphics

which looks a lot better than the before.

Edit: Real coefficients

Here, coeff[n] are the coefficients for the Fourier series in exponential form, but these can be easily converted to the coefficients for the $\cos$ and $\sin$ series, a_n and b_n, by doing something like

a[0] = coeff[0];
a[n_] = Simplify[ComplexExpand[coeff[n] + coeff[-n]]];
b[n_] = Simplify[ComplexExpand[I (coeff[n] - coeff[-n])]];

You can use FourierCoefficient to pre-calculate the Fourier coefficients to arbitrary degree and then use the result very effectively.

There may be some issues with zero-th degree, therefore I excluded this using Piecewise. Here's the main code block:

f[x_] := Piecewise[{
    {-x^3 - 2 x, -2 < x < 0},
    {-1 + x, 0 <= x <= 2}},
    0
];
    Module[{x, fp},
    (* Set parameters so that the integration runs
       from -2 to 2 *)
    fp = {0, -Pi/2};
    fc = Piecewise[{
        {FourierCoefficient[f[x], x, 0, FourierParameters -> fp], # == 0},
        {ComplexExpand@FourierCoefficient[f[x], x, #, FourierParameters -> fp], True}
    }] &;
    fc = Evaluate /@ fc;
];

The output (fc) of this is something very unpleasant to look at; however, there's nothing Fourier-related left there, all the hard math has been done already, and only a bunch of elementary functions remain. fc is now a function of one argument that gives you the $n$-th Fourier coefficient of f in no time.

(* Calculate the first 2001 Fourier coefficients *)
AbsoluteTiming[Table[fc[n], {n, -1000, 1000}]] // First
0.777059 seconds

To convert this back to the function, you have to do the partial sum with your hands, for example here are the second, fourth and eightieth partial sums:

myPartialSums = Table[
    (* The Pi/2 compensates for the custom FourierParameters, see
       documentation of FourierSeries/FourierParameters
       under "more info" *)
    Re[Sum[fc[k] Exp[-Pi/2 I k t], {k, -n, n}]],
    {n, {2, 4, 80}}
];
A very large output has been generated,
but we're luckily not interested in it
anyway but would rather plot it
Plot[
    {f[t]} ~Join~ myPartialSums,
    {t, -10, 10},
    PlotRange -> All, 
    Evaluated -> True, 
    PlotStyle -> {Thick, Automatic, Automatic, Automatic}
]

plot of Fourier series partial sums


If the summation index goes from 1 to n, as you suggested (?) above, and you make the mod to the Plot function I mentioned

s[n_, x_] := 8/4 + 3/(9 π) Sum[(6 (-1)^k)/(k π) Cos[(k π x)/
    2] + (16 (-1)^k + 13)/(π k) Sin[(k π x)/2], {k, 1, n}]

partialsums = Table[s[n, x], {n, 1, 5}];

f[x_] = Piecewise[{{-x^3 - 2 x, -2 < x < 0}, {-1 + x, 0 <= x <= 2}}]

Plot[Evaluate@Tooltip[Append[partialsums, f[x]]], {x, -4 Pi, 4 Pi}]

Then you get the following figure:

enter image description here