echo off;
% NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
% To accompany the text:
% NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
% Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
% This free software is complements of the author.

% Algorithm 6.1 (Differentiation Using Limits).
% Section 6.1, Approximating the Derivative, Page 326
echo on; clc; format long; hold off; clear

% This program finds a numerical approximation for f'(x0).

% The method employed is the limit process.

% The function f(x) is stored in the M-file f.m.
% function z = f(x)
% z = cos(x);

delete f.m
diary f.m; disp('function z = f(x)');...
disp('z = cos(x);');...
diary off;

f(0); % Remark. f.m and difflim.m are used for Algorithm 6.1

pause % Press any key to continue.

clc;
% Place abscissa for differentiation in x0.

% Place the endpoints interval [a,b] about x0 in a and b.

x0 = 0.8;

a = 0;

b = pi/2;

y0 = f(x0);

pause % Press any key to graph the function.

clc; clg;
h = (b-a)/100;
Xs = a:h:b;
Ys = f(Xs);
a = 0;
b = 1.6;
c = 0;
d = 1.1;
axis([a b c d]);...
plot(x0,y0,'or',Xs,Ys,'g');...
hold on;...
plot([a b],[0 0],'b',[0 0],[c d],'b');...
xlabel('x');...
ylabel('y');...
title('y = f(x) and the given point');...
grid;...
axis;...
hold off;...
shg; pause % Press any key to continue.

clc;
% Place abscissa for differentiation in x0.

% Place the tolerance in toler.

x0 = 0.8;

toler = 1e-12;

[H,D,E,n] = difflim('f',x0,toler);

pause % Press any key to see the approximate derivative.

X1 = [a b];
C = [D(n) f(x0)];
Y1 = polyval(C,X1-x0);
axis([a b c d]);...
plot(x0,y0,'or',Xs,Ys,'-g',X1,Y1,'-r');...
hold on;...
plot([-0.05 1.6],[0 0],'b',[0 0],[-0.05 1.05],'b');...
xlabel('x');...
ylabel('y');...
title('The tangent line has slope m = f`(x0).');...
grid;...
axis;...
hold off;...
shg; pause % Press any key to continue.

points = [H;D];
Mx0 = 'Approximating the derivative by finding a limit.';
Mx1 = 'Numerical approximations for ';
Mx2 = ['f`(',num2str(x0),')'];
Mx3 = [Mx1,Mx2];
Mx4 = 'The approximate value of ';
Mx5 = [Mx4,Mx2,' = '];
Mx6 = 'An estimate for the error bound is:';
Mx7 = ' approx. derivative +- error bound';
clc,echo off,diary output,...
disp(''),disp(Mx0),disp(''),disp(Mx3),...
disp(' h D(k)'),disp(points'),...
disp(''),disp(Mx5),disp(D(n)),...
disp(Mx6),disp(Mx7),disp([D(n) E(n)]),diary off,echo on