The Mathematica notebook fourier_series.nb has examples
of such plots.
It's efficient to define a function to make your plot if you've
defined the coefficients as above (also define b0 and the bn
coefficients) :
fs[x_, nmax_] :=
b0 + Sum[an[n]*Sin[n*2 Pi x] + bn[n]*Cos[n*2 Pi x], {n, 1, nmax}]
Plot[{f[x], fs[x, 1], fs[x, 2], fs[x, 3]}, {x, 0, 1}]
For similarities and differences, observe which
expansion does better in different regions (e.g., at the origin),
whether they exhibit the Gibbs overshoot, what values they have
at a discontinuity, whether there is a constant term or not,
and anything else you can think of!
Try extending the range of your plots to -2 to 2. This shows
you the periodic extension of each function.