function [xvect,xdif,fx,nit]=mynewton(x0,nmax,fun,dfun,toll);
% MYNEWTON   Do the newton iteration to find the zeros of the given
% inline scalar function. 
%    [XVEC,XDIF,FX,NIT] = MYNEWTON(X0,NMAX,FUN,DFUN,TOLL)  
%       Input:  x0: starting value
%             nmax: maximum number of iteration
%             toll: tolerance, default is 1e-10
%              fun: given inline function
%             dfun: derivative of the function given as inline
%             function 
%
%    Output: xvect: stores values in all iterations (arg)
%             xdif: difference between two consequent values (arg) 
%               fx: stores function values  in all iterations 
%              nit: number of iteration required
%
%      Examples:
%         The first examples converge to unique zero
%            fun=inline('x^3+x^2-4');  x0=0.3; nmax=100;
%           dfun=inline('3*x^2+2*x');
%
%            fun=inline('x^4+x^3-5*x-12'); 
%           dfun=inline('4*x^3+3*x^2-5');
%
%            fun=inline('x^5+10*x^2-9*x+10');
%           dfun=inline('5*x^4+20*x-9');
%
%         Examples where different starting values converge to
%         different zeros
%            fun=inline('x^2-4*sin(x)');
%           dfun=inline('2*x-4*cos(x)');
%
%            fun=inline('x^3+x^2*(cos(x))^2-4');
%           dfun=inline('(3*x^2+2*x*cos(x)^2-2*x^2*cos(x)*sin(x))');
	
% Author: Bishnu Lamichhane, University of Stuttgart

if (nargin==4)  toll=1e-10; end

err=toll+1; nit=0; xvect=x0;x=x0;fx=feval(fun,x);xdif=[];
  while(nit<nmax & err>toll)
    nit=nit+1; x=xvect(nit);dfx=feval(dfun,x);
    if (dfx==0), err=toll*1e-10;
      disp('Stop for vanishing dfun');
    else,
      xn=x-fx(nit)/dfx;err=abs(xn-x);xdif=[xdif;err];
      x=xn;xvect=[xvect;x];fx=[fx;feval(fun,x)];
    end
  end
n=1:nit;
plot(n, xdif, '-*');
title(['Plot of error with respect to iteration, f(x)=',char(fun)]);
xc = get(gca,'XLim');
yc = get(gca,'YLim');
xc = (xc(1)+xc(2))/2;
text(xc,yc(2)*0.9,['x_{zero} = ',num2str(x),'; x_{start} =' , ...
									 num2str(x0)], 'HorizontalAlignment', 'center'); 
