% 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