function lyapunov(filename)% load in data filedata = load(filename);% calculate number of data pointsN = length(data);% create another array for return mapnextdata = zeros(N,1);for (i=1:N-1)  nextdata(i)=data(i+1);endnextdata(N)=nextdata(N-1);% now plot it (except last point)plot(data(1:N-1),nextdata(1:N-1),'.');% give it a titletitle('return map');display('press any key to continue');% wait for user input ...pause;% calculate the fast fourier transformfftcoeff = fft(data,N);t = 2:N;% plot the absolute value fourier coefficientsplot(t,abs(fftcoeff(2:N)));% give it a titletitle('Fourier Coefficients');display('press any key to continue');% wait for user inputpause;N2 = floor(N/2);N4 = floor(N/4);% find mid point of orbit sequencek=N2;% create space for exponentsexponent = zeros(N4,1);% look at 1/4 of the pointsfor (j=1:N4)    % set distance initially    d = abs(data(k+1)-data(k));    index = k+1;    for (i=2:N-1)        % see if there is a closer point        if (i ~= k) && (abs(data(i)-data(k)))<d            % if so, store index and distance            d = abs(data(i)-data(k));            index = i;        end    end    % write log of quotient as difference of logs to get better accuracy!    if (data(k) ~= data(index)) && (data(k+1) ~= data(index+1))        exponent(j) = log( abs(data(k+1)-data(index+1)))-log(abs(data(k)-data(index)));    end    % repeat with the next point    k = k+1;end% now plot the lyapunov exponentst = 1:N4;lyapunov = exponent(1:N4);exp_avg = 0.0;% find the average value for lyapunov exponentfor (i=1:N4)    exp_avg = exp_avg + exponent(i);end% plot the exponents, the average, and the baselineexp_avg = exp_avg/N4;plot(t,lyapunov,t,0,t,exp_avg);display 'average value for lyapunov exponent is:'exp_avg