%% Chapter 4
% _This is a MATLAB Live Editor version of the *Fundamentals of Numerical Computation* 
% Jupyter Notebooks found at_ <https://tobydriscoll.net/project/fnc/ _https://tobydriscoll.net/project/fnc/_>_._ 
% Example 4.1.1 Bessel Zeros
% In the theory of vibrations of a circular drum, the displacement of the drumhead 
% can be expressed in terms of pure harmonic modes,
% 
% $$J_m(\omega_{k,m} r) \cos(m\theta) \cos(c \omega_{k,m} t),$$
% 
% where $(r,\theta)$ are polar coordinates, $0\le r\le 1$, $t$ is time, $m$ 
% is a positive integer, $c$ is a material parameter, and $J_m$ is a _Bessel function 
% of the first kind_. The quantity $\omega_{k,m}$ is a resonant frequency and 
% is a positive root of the equation
% 
% $$J_m(\omega_{k,m}) = 0,$$
% 
% which states that the drumhead is clamped around the rim. Tabulating approximations 
% to the zeros of Bessel functions has occupied countless mathematician-hours 
% throughout the centuries. 
% 
% Too bad they didn't have MATLAB back in the day. Bessel functions are built 
% in.clear

clf
J3 = @(x) besselj(3,x);
fplot(J3,[0 20]), grid on
xlabel('x'), ylabel('J_3(x)')    % ignore this line
title('Bessel function')     % ignore this line
%% 
% From the graph we see roots near 6, 10, 13, 16, and 19. We use |fzero| to 
% find these roots accurately.

omega = [];
for guess = [6,10,13,16,19]
    omega = [omega;fzero(J3,guess)];
end
omega
%%
hold on, plot(omega,J3(omega),'.')
title('Bessel roots')    
hold off
%% 
% *Example 4.1.2 Rootfinding Condition Number*
% 
% Consider first the function

clear
f  = @(x) (x-1).*(x-2);
%% 
% At the root $r=1$, we have $f'(r)=-1$. If the values of $f$ were perturbed 
% at any point by noise of size, say, $0.05$, we can imagine finding the root 
% of the function as though drawn with a thick line, whose edges we show here.

interval = [0.8 1.2];
fplot(f,interval), grid on, hold on
fplot(@(x) f(x)+0.02,interval,'k')
fplot(@(x) f(x)-0.02,interval,'k')
axis equal                       
xlabel('x'), ylabel('f(x)')      
title('Well-conditioned root')   
%% 
% The possible values for a perturbed root all lie within the interval where 
% the black lines intersect the $x$ axis. The width of that zone is about the 
% same as the vertical distance between the lines.
%% 
% By contrast, consider the function

f  = @(x) (x-1).*(x-1.01);
%% 
% Now $f'(1)=-0.01$, and the graph of $f$ will be much shallower near $x=1$. 
% Look at the effect this has on our thick rendering:

axis(axis), cla
fplot(f,interval)
fplot(@(x) f(x)+0.02,interval,'k')
fplot(@(x) f(x)-0.02,interval,'k')
title('Poorly conditioned root')    
%% 
% The vertical displacements in this picture are exactly the same as before. 
% But the potential _horizontal_ displacement of the root is much wider. In fact, 
% if we perturb the function upward by the amount drawn here, the root disappears 
% entirely!
%% 
% *Example 4.2.1 Fixed Point Sprial*
% 
% Let's convert the roots of a quadratic polynomial $f(x)$ to a fixed point 
% problem.

f = @(x) x.^2 - 4*x + 3.5;
r = roots([1 -4 3.5])
%% 
% We'll define $g(x)=x-f(x)$. Intersections of its graph with the line $y=x$ 
% are fixed points of $g$ and thus roots of $f$. (Only one is shown in the chosen 
% plot range.) 

figure
g = @(x) x - f(x);
fplot(g,[2 3])
hold on, 
plot([2 3],[2 3],'k')
title('Finding a fixed point'), %axis equal   
xlabel('x'), ylabel('y')    
hold off
%% 
% If we evalaute $g(2.1)$, we get a value of almost 2.6. 

x = 2.1;  y = g(x)
%% 
% So $g(x)$ is considerably closer to a fixed point than $x$ was. The value 
% $y=g(x)$ ought to become our new $x$ value! Changing the $x$ coordinate in this 
% way is the same as following a horizontal line over to the graph of $y=x$.

hold on
plot([x y],[y y],'-')
x = y;
hold off
%% 
% Now we can compute a new value for $y$. We leave $x$ alone here, so we travel 
% along a vertical line to the graph of $g$.

hold on
y = g(x)
plot([x x],[x y],'-')
hold off
%% 
% You see that we are in a position to repeat these steps as often as we like. 
% Let's apply them a few times and see the result. 

hold on
for k = 1:5
    plot([x y],[y y],'-'),  x = y;    % y --> new x
    y = g(x);  plot([x x],[x y],'-')  % g(x) --> new y
end
hold off
%% 
% The process spirals in beautifully toward the fixed point we seek. Our last 
% estimate has almost 4 accurate digits.

abs(y-r(1))/r(1)
%% 
% Now let's try to find the other fixed point in the same way. We'll use 1.3 
% as a starting approximation.

figure
fplot(g,[1 2])
hold on, 
plot([1 2],[1 2],'k'), ylim([1 2])
x = 1.3; y = g(x);
for k = 1:5
    plot([x y],[y y],'-'),  x = y;    % y --> new x
    y = g(x);  plot([x x],[x y],'-')  % g(x) --> new y
end
title('No convergence')    
hold off
%% 
% We started near the true fixed point, which is at about 1.293. But this time, 
% the iteration is pushing us _away_ from the correct answer. 
%% 
% *Example 4.2.3 Fixed Point Convergence*

f = @(x) x.^2 - 4*x + 3.5;
r = roots([1 -4 3.5]);
%% 
% Here is the fixed point iteration. This time we keep track of the whole sequence 
% of approximations. 

g = @(x) x - f(x);
x = 2.1; 
for k = 1:12
    x(k+1) = g(x(k));
end
%% 
% It's easiest to construct and plot the sequence of errors. 

figure
err = abs(x-r(1));
semilogy(err,'.-'), axis tight
xlabel('iteration'),  ylabel('error')
title('Convergence of fixed point iteration')
%% 
% It's quite clear that the convergence quickly settles into a linear rate. 
% We could estimate this rate by doing a least-squares fit to a straight line 
% using MATLAB's built-in |polyfit|. Keep in mind that the values for small $k$ 
% should be left out of the computation, as they don't represent the linear trend. 

p = polyfit(5:12,log(err(5:12)),1)
%% 
% The first value is the slope, which must be negative. We can exponentiate 
% it to get the convergence constant $\sigma$.

sigma = exp(p(1))
%% 
% The numerical values of the error should decrease by a factor of $\sigma$ 
% at each iteration. We can check this easily with an elementwise division.

err(9:12) ./ err(8:11)
%% 
% The methods for finding $\sigma$ agree well.
%% 
% *Example 4.3.1 Newton's Method*
% 
% Suppose we want to find a root of the function 

f = @(x) x.*exp(x) - 2;
fplot(f,[0 1.5])
xlabel('x'), ylabel('y')    
set(gca,'ygrid','on')  
title('Objective function')    
%% 
% From the graph, it is clear that there is a root near $x=1$. So we call that 
% our initial guess, $x_1$.

x1 = 1;
f1 = f(x1)
hold on, 
plot(x1,f1,'.','MarkerSize',20)
hold off
%% 
% Next, we can compute the tangent line at the point $\bigl(x_1,f(x_1)\bigr)$, 
% using the derivative.

dfdx = @(x) exp(x).*(x+1);
slope1 = dfdx(x1);
tangent1 = @(x) f1 + slope1*(x-x1);

hold on
fplot(tangent1,[0 1.5],'--')
title('Function and tangent line')
hold off
%% 
% In lieu of finding the root of $f$ itself, we settle for finding the root 
% of the tangent line approximation, which is trivial. Call this $x_2$, our next 
% approximation to the root.

x2 = x1 - f1/slope1

hold on
plot(x2,0,'.','MarkerSize',20)
f2 = f(x2)
title('Root of the tangent')
hold off
%% 
% The residual (value of $f$) is smaller than before, but not zero. So we repeat 
% the process with a new tangent line based on the latest point on the curve.

hold on
fplot(f,[0.8 0.9])
plot(x2,f2,'.','MarkerSize',20)
slope2 = dfdx(x2);
tangent2 = @(x) f2 + slope2*(x-x2);
fplot(tangent2,[0.8 0.9],'--')
x3 = x2 - f2/slope2;
plot(x3,0,'.','MarkerSize',20)
f3 = f(x3)
title('Next iteration')  
hold off
%% 
% We appear to be getting closer to the true root each time. 
%% 
% *Example 4.3.2 Convergence*
% 
% We again look at finding a solution of $xe^x=2$ near $x=1$. To apply Newton's 
% method, we need to calculate values of both the residual function $f$ and its 
% derivative. 

f = @(x) x.*exp(x) - 2;
dfdx = @(x) exp(x).*(x+1);
%% 
% We don't know the exact root, so we use the built-in |fzero| to determine 
% the ``true'' value.

format long,  r = fzero(f,1)
%% 
% We use $x_1=1$ as a starting guess and apply the iteration in a loop, storing 
% the sequence of iterates in a vector.

x = 1;
for k = 1:6
    x(k+1) = x(k) - f(x(k)) / dfdx(x(k));
end
%% 
% Here is the sequence of errors (leaving out the last value, now serving as 
% our exact root). 

format short e
err = x' - r
%% 
% Glancing at the exponents of the errors, they roughly form a neat doubling 
% sequence until the error is comparable to machine precision. We can see this 
% more precisely by taking logs.

format short
logerr = log(abs(err))
%% 
% Quadratic convergence isn't as graphically distinctive as linear convergence.

semilogy(abs(err),'.-')
xlabel('k'), ylabel('|x_k-r|')    
title('Quadratic convergence')    
%% 
% This looks faster than linear convergence, but it's not easy to say more. 
% If we could use infinite precision, the curve would continue to steepen forever.
%% 
% *Example 4.3.3 Usage*
% 
% Suppose we want to solve $e^x=x+c$ for multiple values of $c$. We can create 
% functions for $f$ and $f'$ in each case.

for c = [2 4 7.5 11]
    f = @(x) exp(x) - x - c;
    dfdx = @(x) exp(x) - 1;
    x = newton(f,dfdx,1);  r = x(end);
    fprintf('root with c = %4.1f is %.14f\n',c,r)
end
%% 
% The definition of |f| locks in whatever value is defined for |c| at the moment 
% of definition. Even if we later change the value assigned to |c|, the function 
% is unaffected.

f(r)
c = 100; f(r)
%% 
% This can get a little tricky, because the function is not executed or checked 
% at definition time. You may discover an error using such a definition only later 
% in the code.

clear c
f = @(x) exp(x) - x - c;  % executes OK
c = 1;                    % does not change f
%% 
% A call such as |f(1)| would create an error, since the assignment of |c| did 
% not come until after |f| was defined. 
%% 
% *Example 4.4.1 Secant Method*
% 
% We return to finding a root of the equation $xe^x=2$.

f = @(x) x.*exp(x) - 2;
fplot(f,[0.25 1.25])
set(gca,'ygrid','on')  
xlabel('x'), ylabel('y')    
title('Objective function')    
%% 
% From the graph, it's clear that there is a root near $x=1$. To be more precise, 
% there is a root in the interval $[0.5,1]$. So let us take the endpoints of that 
% interval as _two_ initial approximations. 

x1 = 1;    f1 = f(x1);
x2 = 0.5;  f2 = f(x2);
hold on 
plot([x1 x2],[f1 f2],'.','MarkerSize',20)
title('Two initial guesses')   
hold off
%% 
% Instead of constructing the tangent line by evaluating the derivative, we 
% can construct a linear model function by drawing the line between the two points 
% $\bigl(x_1,f(x_1)\bigr)$ and $\bigl(x_2,f(x_2)\bigr)$. This is called a _secant 
% line_.

slope2 = (f2-f1) / (x2-x1);
secant2 = @(x) f2 + slope2*(x-x2);
%% 
% As before, the next value in the iteration is the root of this linear model. 

hold on
fplot(secant2,[0.25 1.25],'k--')
x3 = x2 - f2/slope2;
f3 = f(x3)

plot(x3,0,'.','MarkerSize',15)
title('Next estimate') 
hold off
%% 
% For the next linear model, we use the line through the two most recent points. 
% The next iterate is the root of that secant line, and so on.

slope2 = (f3-f2) / (x3-x2);
x4 = x3 - f3/slope2;
f4 = f(x4)
%% 
% *Example 4.4.2 Convergence*
% 
% We check the convergence of the secant method from the previous example.

f = @(x) x.*exp(x) - 2;
x = secant(f,1,0.5);
%% 
% We don't know the exact root, so we use the built-in |fzero| to get a substitute.

r = fzero(f,1);
%% 
% Here is the sequence of errors. 

format short e
err = r - x(1:end-1)'
%% 
% It's not so easy to see the convergence rate by looking at these numbers. 
% But we can check the ratios of the log of successive errors. 

logerr = log(abs(err));
ratios = logerr(2:end) ./ logerr(1:end-1)
%% 
% It seems to be heading toward a constant ratio of about 1.6 by the time it 
% quits.
%% 
% *Example 4.4.3 Inverse Interpolation*
% 
% Here we look for a root of $x+\cos(10x)$ that is close to 1. 

f = @(x) x + cos(10*x);
interval = [0.5,1.5];
fplot(f,interval)
set(gca,'ygrid','on'), axis(axis)   
title('Objective function')    
xlabel('x'), ylabel('y')    
r = fzero(f,1)
%% 
% We choose three values to get the iteration started.

x = [0.8 1.2 1]';
y = f(x);
hold on
plot(x,y,'.','MarkerSize',15)
title('Three initial points')   
hold off
%% 
% If we were using "forward" interpolation, we would ask for the polynomial 
% interpolant of $y$ as a function of $x$. But that parabola has no real roots. 

c = polyfit(x,y,2);    % coefficients of interpolant
q = @(x) polyval(c,x);
hold on
fplot(q,interval,'--')
fplot(f,interval)
title('Parabola model') 
axis([.5 1.5 -.5 2.5]);
hold off
%% 
% To do inverse interpolation, we swap the roles of $x$ and $y$ in the interpolation.

cla, 
fplot(f,interval), 
hold on
plot(x,y,'.','MarkerSize',15)     
c = polyfit(y,x,2);    % coefficients of interpolating polynomial
q = @(y) polyval(c,y);

fplot(q,@(y) y,ylim,'--')    % plot x=q(y), y=y
title('Sideways parabola') 
fplot(f,interval)
hold off
%% 
% We seek the value of $x$ that makes $y$ zero. This means evaluating $q$ at 
% zero. 

x = [x; q(0)];
y = [y; f(x(end))]
%% 
% We repeat the process a few more times.

for k = 4:8
    c = polyfit(y(k-2:k),x(k-2:k),2);
    x(k+1) = polyval(c,0);
    y(k+1) = f(x(k+1));
end
%% 
% Here is the sequence of errors.

format short e
err = x - r
%% 
% The error seems to be superlinear, but subquadratic.

logerr = log(abs(err));
ratios = logerr(2:end) ./ logerr(1:end-1)
%% 
% *Example 4.5.2 Nonlinear Systems*
% 
% At the end of this file you will find functions that define a system of three 
% equations and its Jacobian matrix. (These functions could be written as separate 
% files, or embedded within another function as we have done here.) Our initial 
% guess at a root is the origin.

x = [0;0;0];   % column vector!
%% 
% We need the value (residual) of the nonlinear function, and its Jacobian, 
% at this value for $\mathbf{x}$. 

f = nlvalue(x);
J = nljac(x);
%% 
% We solve for the Newton step and compute the new estimate. 

s = -(J\f);
x(:,2) = x(:,1) + s;
%% 
% Here is the new residual.

format short e
f(:,2) = nlvalue(x(:,2))
%% 
% We don't seem to be especially close to a root. Let's iterate a few more times. 

for n = 2:7
    f(:,n) = nlvalue(x(:,n));
    s = -( nljac(x(:,n)) \ f(:,n) );
    x(:,n+1) = x(:,n) + s;
end
%% 
% We find the infinity norm of the residuals. 

residual_norm = max(abs(f),[],1)'   % max in dimension 1
%% 
% We don't know an exact answer, so we will take the last computed value as 
% its surrogate. 

r = x(:,end);
x = x(:,1:end-1);
%% 
% The following will subtract r from every column of x.

z = size(x);
e = x - r*ones(1,z(2));
%% 
% Now we take infinity norms and check for quadratic convergence. 

logerrs = log(abs( max(e,[],1)' ));
ratios = logerrs(2:end) ./ logerrs(1:end-1)
%% 
% For a brief time, we see ratios around 2, but then the limitation of double 
% precision makes it impossible for the doubling to continue. 
%% 
% *Example 4.5.3 Usage*
% 
% We repeat a previous example. The system is again defined by its residual 
% and Jacobian, but this time we implement them as a single function (found at 
% the end of the file).
% 
% Our initial guess is the origin. The output has one column per iteration.

x1 = [0;0;0];   % column vector!
x = newtonsys(@nlsystem,x1);
[~,num_iter] = size(x)
%% 
% The last column contains the final Newton estimate. We'll compute the residual 
% there in order to check the quality of the result.

r = x(:,end)
back_err = norm(nlsystem(r))
%% 
% Let's use the convergence to the first component of the root as a proxy for 
% the convergence of the vectors.

log10( abs(x(1,1:end-1)-r(1)) )'
%% 
% The exponents approximately double, as is expected of quadratic convergence. 
%% 
% *Example 4.6.1 Levenberg*
% 
% We repeat a previous example. We need to code only the function defining the 
% system, and not its Jacobian. The function is found at the end of this file.
% 
% In all other respects usage is the same as for the |newtonsys| function. 

x1 = [0;0;0];   
x = levenberg(@nlsystem2,x1);
%% 
% It's always a good idea to check the accuracy of the root, by measuring the 
% residual (backward error). 

r = x(:,end)
backward_err = norm(nlsystem2(r))
%% 
% Looking at the convergence of the first component, we find a subquadratic 
% convergence rate, just as with the secant method.

log10( abs(x(1,1:end-1)-r(1)) )'
%% 
% *Example 4.7.1 Michaelis--Menten Fit*
% 
% Inhibited enzyme reactions often follow what are known as *Michaelis--Menten* 
% kinetics, in which a reaction rate $v$ follows a law of the form
% 
% $$v(x) = \frac{V x}{K_m + x},$$
% 
% where $x$ is the concentration of a substrate. The real values $V$ and $K_m$ 
% are parameters that are free to fit to data. For this example we cook up some 
% artificial data with $V=2$ and $K_m=1/2$.

clear
global x_data y_data   % FNC did not declare global variables

m = 25;
x_data = linspace(0.05,6,m)';
y_data = 2*x_data./(0.5+x_data);                   % exactly on the curve
y_data = y_data + 0.15*cos(2*exp(x_data/16).*x_data);   % noise added
plot(x_data,y_data,'.','MarkerSize',20)
title('Data')   
xlabel('x'), ylabel('v') 
%% 
% The idea is to pretend that we know nothing of the origins of this data and 
% use nonlinear least squares on the misfit to find the parameters in the theoretical 
% model function $v(x)$. Note in the Jacobian that the derivatives are _not_ with 
% respect to $x$, but with respect to the two parameters, which are contained 
% in the vector |c|. The function coding the equations and Jacobian is found at 
% the end of this file. 

c1 = [1; 0.75];
c = newtonsys(@misfit,c1);
V = c(1,end),  Km = c(2,end)   % final values
model = @(x) V*x./(Km+x);
%% 
% The final values are close to the noise-free values of $V=2$, $K_m=0.5$ that 
% we used to generate the data. We can calculate the amount of misfit at the end, 
% although it's not completely clear what a "good" value would be. Graphically, 
% the model looks reasonable.

final_misfit_norm = norm(model(x_data)-y_data) 
hold on
fplot( model, [0 6] )
title('Michaelis-Menten fitting')   
xlabel('concentration'), ylabel('reaction rate')    
hold off
%% 
% For this model, we also have the option of linearizing the fit process. Rewrite 
% the model as $1/v= (a/x)+b$ for the new parameters $\alpha=K_m/V$ and $\beta=1/V$. 
% This corresponds to the misfit function whose entries are
% 
% $$f_i(\alpha,\beta) = \alpha \cdot \frac{1}{x_i} + \beta - \frac{1}{y_i},$$
% 
% for $i=1,\ldots,m$. Although the misfit is nonlinear in $x$ and $y$, it's 
% linear in the unknown parameters $\alpha$ and $\beta$, and so can be posed and 
% solved as a linear least-squares problem.

A = [ x_data.^(-1), x_data.^0 ];  u = 1./y_data;
z =  A\u;
alpha = z(1);  beta = z(2);
%% 
% The two fits are different, because they do not optimize the same quantities. 

linmodel = @(x) 1 ./ (beta + alpha./x);
final_misfit_linearized = norm(linmodel(x_data)-y_data)
hold on
fplot(linmodel, [0 6] )
hold off
%% 
% _*Below are function defined for the above examples. All functions need to 
% be at the end of a Live Editor script.*_

function [f,J] = misfit(c)

    global x_data y_data
    V = c(1);   Km = c(2);
    f = V*x_data./(Km+x_data) - y_data;
    J(:,1) = x_data./(Km+x_data);            % d/d(V)
    J(:,2) = -V*x_data./(Km+x_data).^2;      % d/d(Km)
end

function x = newton(f,dfdx,x1)
% NEWTON   Newton's method for a scalar equation.
% Input:
%   f        objective function 
%   dfdx     derivative function
%   x1       initial root approximation
% Output       
%   x        vector of root approximations (last one is best)

% Operating parameters.
funtol = 100*eps;  xtol = 100*eps;  maxiter = 40;

x = x1;  
y = f(x1);
dx = Inf;   % for initial pass below
k = 1;

while (abs(dx) > xtol) && (abs(y) > funtol) && (k < maxiter)
  dydx = dfdx(x(k));
  dx = -y/dydx;           % Newton step
  x(k+1) = x(k) + dx;

  k = k+1;
  y = f(x(k));
end

if k==maxiter, warning('Maximum number of iterations reached.'), 
end
end

function x = secant(f,x1,x2)
% SECANT   Secant method for a scalar equation.
% Input:
%   f        objective function 
%   x1,x2    initial root approximations
% Output       
%   x        vector of root approximations (last is best)

% Operating parameters.
funtol = 100*eps;  xtol = 100*eps;  maxiter = 40;

x = [x1 x2];
dx = Inf;  y1 = f(x1);
k = 2;  y2 = 100;

while (abs(dx) > xtol) && (abs(y2) > funtol) && (k < maxiter)
  y2 = f(x(k));
  dx = -y2 * (x(k)-x(k-1)) / (y2-y1);   % secant step
  x(k+1) = x(k) + dx;
  
  k = k+1;
  y1 = y2;    % current f-value becomes the old one next time    
end

if k==maxiter, warning('Maximum number of iterations reached.'), 
end
end

function x = newtonsys(f,x1)
% NEWTONSYS   Newton's method for a system of equations.
% Input:
%   f        function that computes residual and Jacobian matrix
%   x1       initial root approximation (n-vector)
% Output       
%   x        array of approximations (one per column, last is best)

% Operating parameters.
funtol = 1000*eps;  xtol = 1000*eps;  maxiter = 40;

x = x1(:);  
[y,J] = f(x1);
dx = Inf;
k = 1;

while (norm(dx) > xtol) && (norm(y) > funtol) && (k < maxiter)
  dx = -(J\y);   % Newton step
  x(:,k+1) = x(:,k) + dx;

  k = k+1;
  [y,J] = f(x(:,k));
end

if k==maxiter, warning('Maximum number of iterations reached.')
end
end

%% 
% 

function f = nlsystem2(x)
    f = zeros(3,1);   % ensure a column vector output
    f(1) = exp(x(2)-x(1)) - 2;
    f(2) = x(1)*x(2) + x(3);
    f(3) = x(2)*x(3) + x(1)^2 - x(2);
end

function [f,J] = nlsystem(x)
    f = zeros(3,1);   % ensure a column vector output
    f(1) = exp(x(2)-x(1)) - 2;
    f(2) = x(1)*x(2) + x(3);
    f(3) = x(2)*x(3) + x(1)^2 - x(2);
    J(1,:) = [-exp(x(2)-x(1)),exp(x(2)-x(1)), 0];
    J(2,:) = [x(2), x(1), 1];
    J(3,:) = [2*x(1), x(3)-1, x(2)];
end
%% 
% This function defines the system of nonlinear equations.

function f = nlvalue(x)
    f = zeros(3,1);   % ensure a column vector output
    f(1) = exp(x(2)-x(1)) - 2;
    f(2) = x(1)*x(2) + x(3);
    f(3) = x(2)*x(3) + x(1)^2 - x(2);
end
%% 
% Here is a function that computes the Jacobian matrix.

function J = nljac(x)
    J(1,:) = [-exp(x(2)-x(1)),exp(x(2)-x(1)), 0];
    J(2,:) = [x(2), x(1), 1];
    J(3,:) = [2*x(1), x(3)-1, x(2)];
end
%% 
% Quasi-Newton method for nonlinear systems.

function x = levenberg(f,x1,tol)
% LEVENBERG   Quasi-Newton method for nonlinear systems.
% Input:
%   f         objective function 
%   x1        initial root approximation
%   tol       stopping tolerance (default is 1e-12)
% Output       
%   x         array of approximations (one per column)

% Operating parameters.
if nargin < 3, tol = 1e-12; end
ftol = tol;  xtol = tol;  maxiter = 40;

x = x1(:);     fk = f(x1);
k = 1;  s = Inf;        
Ak = fdjac(f,x(:,1),fk);   % start with FD Jacobian
jac_is_new = true;
I = eye(length(x));

lambda = 10; 
while (norm(s) > xtol) && (norm(fk) > ftol) && (k < maxiter)
    % Compute the proposed step.
    B = Ak'*Ak + lambda*I;
    z = Ak'*fk;
    s = -(B\z);

    xnew = x(:,k) + s;   fnew = f(xnew);
    
    % Do we accept the result?
    if norm(fnew) < norm(fk)    % accept
        y = fnew - fk;
        x(:,k+1) = xnew;  fk = fnew;  
        k = k+1;
        
        lambda = lambda/10;  % get closer to Newton
        % Broyden update of the Jacobian.
        Ak = Ak + (y-Ak*s)*(s'/(s'*s));
        jac_is_new = false;
    else                       % don't accept
        % Get closer to steepest descent.
        lambda = lambda*4;
        % Re-initialize the Jacobian if it's out of date.
        if ~jac_is_new
            Ak = fdjac(f,x(:,k),fk);
            jac_is_new = true;
        end
    end
end

if (norm(fk) > 1e-3)
    warning('Iteration did not find a root.')
end

end

function J = fdjac(f,x0,y0)
% FDJAC   Finite-difference approximation of a Jacobian.
% Input:
%   f        function to be differentiated
%   x0       evaluation point (n-vector)
%   y0       value of f at x0 (m-vector)
% Output       
%   J        approximate Jacobian (m-by-n)

delta = sqrt(eps);   % FD step size
m = length(y0);  n = length(x0);
J = zeros(m,n);
I = eye(n);
for j = 1:n
    J(:,j) = ( f(x0+delta*I(:,j)) - y0) / delta;
end

end