The "tol" value is the initial amount of perturbation to apply in the first iteration. It is the range of the stock price end points times the multiplier (in this case 2/3). Play with the multiplier to get more or less perturbation.

The reason for the exponential is to reduce the magnitude for each successive iteration. I suppose a different function could be used but an exponential reduction replicates the fractal nature (self-similarity upon zooming in).

Thanks for the submission. I had a different way to do this... maybe it will help? If not, no problem. Feel free to incorporate any of my code if you like.

Here is mine (only works on vectors):

function dy = diff2(y)
% DIFF2 computes the first derivative while preserving the number of
% elements. dx is considered to be 1. So, divide this result by your grid
% spacing to get the correct value of the derivative.
%
% DIFF2(y) returns the central difference of y. y is a vector of at least
% three elements.
%
% This scheme is second-order accurate (central differences).
%
% EXAMPLE:
% dx = 0.001; x = 0:dx:2*pi;
% y = sin(x); dydx = diff2(y)/dx;
% plot(x,y,x,dydx);
%
% Created by Mike Hanchak, October 2010
%
% See also DIFF

% Check inputs.
[m,n] = size(y);
if m*n < 3 || m*n > max([m,n]) || nargin ~= 1
error('Input must be a vector of at least three elements.')
end

% Force data into a row vector.
y = y(:)';

% Calculate central difference for interior terms.
temp = (y(3:end)-y(1:end-2))/2;

% Concatenate derivatives at endpoints using second order forward and
% backward differences.
dy = [ (-3*y(1)+4*y(2)-y(3))/2, temp, (3*y(end)-4*y(end-1)+y(end-2))/2];