
%% Example From Global Behavior with eqn no. (1.37)

%% Preparation

% define function and its derivative
f=@(x) x.*(x.^2-1).*(x-3).*exp(-0.5*(x-1).^2);
df=@(x) -exp(-(1/2)*(x-1).^2).*(x.^5-4*x.^4-2*x.^3+13*x.^2-x-3);

% preliminary plotting
x=-3:0.01:5;
plot(x,f(x), '-b', x, 0*x, '-k')

% maximum interation number
maxIt=500;

% define initial values to be examined
nrInitVals=5000;
x0Vals=linspace(-3,5,nrInitVals);

% where to store data
iterates=zeros(maxIt,nrInitVals);

%% Iterate
for j=1:nrInitVals
    
    xk=x0Vals(j);
    for k=1:maxIt

        xk_p1=xk-f(xk)/df(xk);
        iterates(k,j)=xk_p1;

        xk=xk_p1;
        
    end    
end    

%% Visualize

% extract final values (that are not NaNs)
finalVals=zeros(size(x0Vals));
for j=1:nrInitVals
    
    jIt=iterates(:,j);
    jItReal=jIt(isnan(jIt)==0);
    
    finalVals(j)=jItReal(end);
    
end 

% to identify -inf and inf I call them "-3" and "5"
finalVals(finalVals<-2)=-3;
finalVals(finalVals>4)=5;

% all limits (I checked that they really only consist of the zeros and "\pm \inf=-3 and 5")
endVals=union(finalVals,[]);

% set colors to be used
colorVals={'k','r','b','g','m','k'};

% plot the behavior

clf; hold on
for i=1:length(endVals)
    
    iFinal=endVals(i);
    which=finalVals==iFinal;
    
    X=[x0Vals(which);finalVals(which)];
    Y=[ones(size(x0Vals(which))); (-2)*ones(size(x0Vals(which)))];
    
    plot(X,Y, colorVals{i}, 'LineWidth',0.1)
end    

ax=gca;
ax.YTick=[-2; 1];
ax.YTickLabels={'Limit',  'Start Value'};
ax.XTick=[-3 -1 0 1 3 5]';
ax.XTickLabels={'-Inf',  '-1', '0', '1',  '3',  'Inf'};
ax.XGrid='on';
ax.GridAlpha=1;
ax.Layer='top';
ax.FontSize=14;
title('Global Behavior of Newton Method')

