% Solve the nonlinear equation
%        x * (x^2-1) * (x-3) * exp(-0.5(x-1)^2) = 0
% using Newton's method.  This is equation has roots at x=-1, x=0,
% x=1, x=3. The solution Newton's method finds depends on the initialization.

clear all;
clc;
xx = -2:0.01:4;
ff = @(xx)(xx.*(xx.^2-1).*(xx-3).*exp(-0.5*(xx-1).^2));
ffprime = @(xx)((xx.^2-1).*(xx-3).*exp(-0.5*(xx-1).^2) + ...
    2*xx.*xx.*(xx-3).*exp(-0.5*(xx-1).^2) + ...
    xx.*(xx.^2-1).*exp(-0.5*(xx-1).^2) - ...
    xx.*(xx.^2-1).*(xx-3).*exp(-0.5*(xx-1).^2).*(-xx-1));
plot(xx,ff(xx), '-k', 'Linewidth', 2); hold on;
plot(xx,0.*xx, '--r', 'Linewidth', 2);
plot(xx,ffprime(xx), '--g', 'Linewidth', 2);
legend('f(x)', 'y=0', 'fprime');
axis([-2,4,-3,1.5]);

% initialization
x(1) = 3.5;

fprintf('\nNewtons method \n');
fprintf('%10s    %12s   \n', 'Iteration', 'Value');
fprintf('%10d    %.16f  \n', 1, x(1));
for i=1:100
    x(i+1) = x(i) - ff(x(i))/ffprime(x(i));
    fprintf('%10d    %.16f\n', i+1, x(i+1));

    if (abs(ff(x(i+1)))<1e-14)
	fprintf('Converged after %d iterations.\n', i)
	break;
    end
end