% Interval endpoints
c = 0;
d = 1;

% Interpolation points
x = [1,0,-1,0];
n = length(x);   assert(mod(n,2) == 0);
t = c+(d-c)*(0:n-1)/n;

% Compute expansion coefficitnes (FFT)
y = fft(x)/sqrt(n);
a = real(y);
b = imag(y);

% Plot the interpolation function
P_n = @(tj) (a(1) + 2*sum( ...
    a(2:n/2).*cos(2*pi* (tj-c)/(d-c) *(1:n/2-1) ) ...
    - ...
    b(2:n/2).*sin(2*pi* (tj-c)/(d-c) *(1:n/2-1) ) ...
    ) ...
    + a(n/2+1)*cos(2*pi* (tj-c)/(d-c) *(n/2) ) ...
    - b(n/2+1)*sin(2*pi* (tj-c)/(d-c) *(n/2) ) )/sqrt(n);

N = 100;
that = linspace(c,d,N);
Phat = arrayfun(P_n,that);

plot(t,x,'o',that,Phat,'-',...
     'LineWidth',3,'MarkerSize',10);          % plot data points and interpolant
xlabel('t');
ylabel('P');