function df = calc_Jac(x,Param)
% function to calculate the Jacobian matrix for a set of nonlinear
% algebraic equations.  A simple central differencing scheme is used.

% INPUT
% x         a vector of independent variables
% Param     data structure containing parameters

% OUTPUT
% df    the Jacobian matrix

% ASSIGNMENT
diff = 0.01;  % the default differencing percentage for each variable

% create the Jacobian based on the size of input vector x
n = length(x);
df = zeros(n,n);

% if nargin == 1  % set default differencing at 1% for each variable
    xdif = zeros(n,1)+ diff;
% end

% loop to calculate successive columns of the Jacobian
for icols = 1:n
    x_plus = x;
    x_minus = x;
    % create a difference at element icols in the x vector
    x_plus(icols)=x(icols)*(1 + xdif(icols));
    x_minus(icols) = x(icols)*(1 - xdif(icols));
    df(:,icols) = (calc_f(x_plus,Param) - calc_f(x_minus,Param))...
          /(2*xdif(icols)*x(icols));
end

% Dr. Barry S. Johnston, Copyright 2004.