%% Chapter 4
% This is a MATLAB Live Editor version of the Fundamentals of Numerical Computation Jupyter Notebooks found at 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