%% Mfile name
%        mydiff.m

%% Revised: 
%       October 2, 2007

%% Author
%       Autar Kaw
%       Section: All
%       Semester: Fall 2007

%% Purpose
%       Given the function f, find the value of the first derivative of the
%       function at x within a prespecified percentage tolerance tol and a
%       starting step size of dx.
%       We are using the forward divided difference method 
%       f'(x) = [f(x+dx)-f(x)]/dx
%       to find the approximate value of the derivative.

%% Usage
%       fp=mydiff(f,x,dx,tol)
% Input variables
%       f = function
%       x = point at which derivative of f is sought
%       dx = starting step size
%       tol = prespecified percentage tolerance
% Output variables
%       fp = approximate value of derivative of f at point x 

% Keyword
%       Differentiation
%       Forward Divided Difference
%% Code
function [fp]=mydiff(f,x,dx,tol)
%Initializing the step size to a different variable 
% as input variable should not be placed on the LHS of an assignment
dx1=dx;

%Calculating initial value of the derivative of the function at x
fpold=(f(x+dx1)-f(x))/dx1;

%Initializing the absolute relative approximate error to be greater than 
% the presepcified tolerance
% so that the while loop works the first time
absea=2*tol;

% Keep refining the value of the derivative of the function till tolerance is met
while (absea>tol) 
%Halving the step size after each iteration
         dx1=dx1/2.0;
         fpnew=(f(x+dx1)-f(x))/dx1;
%Calculating the absolute relative approximate error
         absea=abs((fpnew-fpold)/fpnew)*100.0;
% Reinitializing old value with a new value
         fpold=fpnew;
end

%Returning the derivative of the function within the pre-specified tolerance
fp = fpnew;