% Chapter 10
% Comparison of Basic Numerical Differentiation Formulas

clear all
close all
clc

f = @(x) 1./(1+25.*x.^2);   % Runge Function
df = @(x) -50.*x./((1+25.*x.^2).^2);
ddf = @(x) (5000.*x.^2)./((1+25.*x.^2).^3) - 50./((1+25.*x.^2).^2);
dddf = @(x) (-750000.*x.^3)./((1+25.*x.^2).^4) + 15000.*x./((1+25.*x.^2).^3);
ddddf = @(x) (150000000.*x.^4)./((1+25.*x.^2).^5)-(4500000.*x.^2)./((1+25.*x.^2).^4)+(15000)./((1+25.*x.^2).^3);
a = -1;
b = 1;

h = 0.1;

xplot = linspace(a,b);
xp = [a:h:b];
yp = f(xp);

figure
subplot(5,3,2)
hold on
plot(xplot,f(xplot),'-')
plot(xp,f(xp),'ok','MarkerFaceColor','k','MarkerSize',3)
xlim([-1,1])
ylim([0,1])
plot([0 0],ylim,'-k')
plot(xlim,[0 0],'-k')
title('Original Function and Data Points')

%--------------------------------------------------------------------------
% First Derivative FFD
n = 1;
[xffd1,ffd1] = ffd(xp,yp,n);
subplot(5,3,4)
hold on
plot(xplot,df(xplot),'-')
plot(xffd1,ffd1,'ok','MarkerFaceColor','k','MarkerSize',3)
xlim([-1,1])
ylim([-4,4])
plot([0 0],ylim,'-k')
plot(xlim,[0 0],'-k')
title('First Derivative FFD1h1')

% First Derivative BFD
n = 1;
[xbfd1,bfd1] = bfd(xp,yp,n);
subplot(5,3,5)
hold on
plot(xplot,df(xplot),'-')
plot(xbfd1,bfd1,'ok','MarkerFaceColor','k','MarkerSize',3)
xlim([-1,1])
ylim([-4,4])
plot([0 0],ylim,'-k')
plot(xlim,[0 0],'-k')
title('First Derivative BFD1h1')

% First Derivative CFD
n = 1;
[xcfd1,cfd1] = cfd(xp,yp,n);
subplot(5,3,6)
hold on
plot(xplot,df(xplot),'-')
plot(xcfd1,cfd1,'ok','MarkerFaceColor','k','MarkerSize',3)
xlim([-1,1])
ylim([-4,4])
plot([0 0],ylim,'-k')
plot(xlim,[0 0],'-k')
title('First Derivative CFD1h2')

%--------------------------------------------------------------------------
% Second Derivative FFD
n = 2;
[xffd2,ffd2] = ffd(xp,yp,n);
subplot(5,3,7)
hold on
plot(xplot,ddf(xplot),'-')
plot(xffd2,ffd2,'ok','MarkerFaceColor','k','MarkerSize',3)
xlim([-1,1])
ylim([-60,20])
plot([0 0],ylim,'-k')
plot(xlim,[0 0],'-k')
title('Second Derivative FFD2h1')

% Second Derivative BFD
n = 2;
[xbfd2,bfd2] = bfd(xp,yp,n);
subplot(5,3,8)
hold on
plot(xplot,ddf(xplot),'-')
plot(xbfd2,bfd2,'ok','MarkerFaceColor','k','MarkerSize',3)
xlim([-1,1])
ylim([-60,20])
plot([0 0],ylim,'-k')
plot(xlim,[0 0],'-k')
title('Second Derivative BFD2h1')

% Second Derivative CFD
n = 2;
[xcfd2,cfd2] = cfd(xp,yp,n);
subplot(5,3,9)
hold on
plot(xplot,ddf(xplot),'-')
plot(xcfd2,cfd2,'ok','MarkerFaceColor','k','MarkerSize',3)
xlim([-1,1])
ylim([-60,20])
plot([0 0],ylim,'-k')
plot(xlim,[0 0],'-k')
title('Second Derivative CFD2h2')

%--------------------------------------------------------------------------
% Third Derivative FFD
n = 3;
[xffd3,ffd3] = ffd(xp,yp,n);
subplot(5,3,10)
hold on
plot(xplot,dddf(xplot),'-')
plot(xffd3,ffd3,'ok','MarkerFaceColor','k','MarkerSize',3)
xlim([-1,1])
ylim([-600,600])
plot([0 0],ylim,'-k')
plot(xlim,[0 0],'-k')
title('Third Derivative FFD3h1')

% Third Derivative BFD
n = 3;
[xbfd3,bfd3] = bfd(xp,yp,n);
subplot(5,3,11)
hold on
plot(xplot,dddf(xplot),'-')
plot(xbfd3,bfd3,'ok','MarkerFaceColor','k','MarkerSize',3)
xlim([-1,1])
ylim([-600,600])
plot([0 0],ylim,'-k')
plot(xlim,[0 0],'-k')
title('Third Derivative BFD3h1')

% Third Derivative CFD
n = 3;
[xcfd3,cfd3] = cfd(xp,yp,n);
subplot(5,3,12)
hold on
plot(xplot,dddf(xplot),'-')
plot(xcfd3,cfd3,'ok','MarkerFaceColor','k','MarkerSize',3)
xlim([-1,1])
ylim([-600,600])
plot([0 0],ylim,'-k')
plot(xlim,[0 0],'-k')
title('Third Derivative CFD3h2')

%--------------------------------------------------------------------------
% Fourth Derivative FFD
n = 4;
[xffd4,ffd4] = ffd(xp,yp,n);
subplot(5,3,13)
hold on
plot(xplot,ddddf(xplot),'-')
plot(xffd4,ffd4,'ok','MarkerFaceColor','k','MarkerSize',3)
xlim([-1,1])
ylim([-7000,16000])
plot([0 0],ylim,'-k')
plot(xlim,[0 0],'-k')
title('Fourth Derivative FFD4h1')

% Fourth Derivative BFD
n = 4;
[xbfd4,bfd4] = bfd(xp,yp,n);
subplot(5,3,14)
hold on
plot(xplot,ddddf(xplot),'-')
plot(xbfd4,bfd4,'ok','MarkerFaceColor','k','MarkerSize',3)
xlim([-1,1])
ylim([-7000,16000])
plot([0 0],ylim,'-k')
plot(xlim,[0 0],'-k')
title('Fourth Derivative BFD4h1')

% Fourth Derivative CFD
n = 4;
[xcfd4,cfd4] = cfd(xp,yp,n);
subplot(5,3,15)
hold on
plot(xplot,ddddf(xplot),'-')
plot(xcfd4,cfd4,'ok','MarkerFaceColor','k','MarkerSize',3)
xlim([-1,1])
ylim([-7000,16000])
plot([0 0],ylim,'-k')
plot(xlim,[0 0],'-k')
title('Fourth Derivative CFD4h2')