%% Fourier coefficients using symbolic integration.

tic;
f = @ (x) x.^2;
syms x k;

% Calculate coefficients symbolically
zeroth_term = int(f(x),x,-pi,pi)/(2*pi);
a_k = int(f(x)*cos(k.*x),x,-pi,pi) / pi;
b_k = int(f(x)*sin(k.*x),x,-pi,pi) / pi;

x = linspace(-pi,pi);
want_terms = 20;

y = zeros(1, length(x)); % Reserve memory
for n = 1:want_terms
    y = y + subs(a_k, k, n) * cos(n*x) + subs(b_k, k, n) * sin(n*x);
end
y = y + zeroth_term;

plot(x,y,'r',x,f(x),'k','linewidth',2);
toc

%% Fourier coefficients using numerical integration.

tic;
f = @ (x) x.^2;

zeroth_term = integral(f,-pi,pi)/(2*pi);

x = linspace(-pi,pi);
want_terms = 20;

y = zeros(1, length(x));  % Reserve memory, for speed.
for n = 1:want_terms,
    a_n = integral(@(xx) f(xx).*cos(n*xx), -pi, pi) / pi;
    b_n = integral(@(xx) f(xx).*sin(n*xx), -pi, pi) / pi;
    y = y + a_n * cos(n*x) + b_n * sin(n*x);
end
y = y + zeroth_term;

plot(x,y,'r',x,f(x),'k','linewidth',2);
toc