% FD4Demo.m illustrates convergence of fourth-order finite differences
% for derivative estimation of the periodic function y = exp(sin(t))
% Taken from Trefethen: Spectral Methods in Matlab

clear all
clc
close all
warning off
scrnsz = get(0, 'Screensize');

% For various N, set up grid in [-pi,pi] and function y(x):
Nvec = 2.^(3:12);
figure('Position', [5 10 scrnsz(3)-5 scrnsz(4)-60])
for N = Nvec
    h = 2*pi/N; t = -pi + (1:N)'*h;
    y = exp(sin(t)); yprime = cos(t).*y;

    % Construct sparse fourth-order differentiation matrix:
    e = ones(N,1);
    D =   sparse(1:N,[2:N 1],2*e/3,N,N)...
        - sparse(1:N,[3:N 1 2],e/12,N,N);
    D = (D-D')/h;

    % Plot max(abs(D*y-yprime)):
    error = norm(D*y-yprime,inf);
    loglog(N,error,'.','markersize',15), hold on
end
grid on, xlabel N, ylabel error
title('Convergence of fourth-order finite differences')
semilogy(Nvec,Nvec.^(-4),'--') 
text(105,5e-8,'N^{-4}','fontsize',18)
hold off