% 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);

% Create the "hatted" coefficients for the "fine" mesh:
N = 128;
yhat = zeros(1,N);
% Transfer them from y
% low frequencies
yhat(1:n/2)   = sqrt(N/n)*y(1:n/2);
% "middle coefficient"
yhat(n/2+1)   = sqrt(N/n)*y(n/2+1)/2;
yhat(N-n/2+1) = sqrt(N/n)*y(n/2+1)'/2;
% high frequencies
yhat(N-n/2+2:N) = sqrt(N/n)*y(n/2+2:n);

% Plot the interpolation function
that = c+(d-c)*(0:N-1)/N;
Phat = real(ifft(yhat))*sqrt(N);

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