% Halley's iteration for n=1 (cubic convergence)
% K. Ming Leung, 02/13/03

clear all; 
tol=eps;    % set the tolerance
x=3;        % initial guess
dx=1;       % set initial change in x 
            % (need this since Matlab have no do-while loop)
result=[];  % initalize an empty matrix to store intermediate results
loop=0;     % loop count
while abs(dx)> tol;     % may use better termination condition
    fx=x*x-4*sin(x);    % the function of interest
    fp=2*x-4*cos(x);    % its first derivative
    fpp=2+4*sin(x);     % its second derivative
    dx=-2*(fx*fp)/(2*fp*fp-fx*fpp); % change in x
    result=[result; [loop x fx fp dx]]; % store intermediate results
    x=x+dx;             % the next value of x
    loop=loop+1;        % increase loop counter
end;
display(result);        % display final result

