Simple modification of Newton-Raphson method to Broyden's method? Faster convergence?

Hello everybody, first and foremost thank you for your time. I am attempting to speed up my program, which solves sets of non-linear equations (mass and energy balances within a chemical reactor). I am currently using the Newton-Raphson method, which works well but is quite slow (I am using a waiting factor to make the solver more stable).

My question is, can someone offer some guidance on how to modify my existing code to utilize Broyden's method, which as I understand it, modifies the Jacobian matrix in relation to the change in f(x) for each iteration, rather than calculating a new Jacobian each time. Any input, in relatively simple terms would be greatly appreciated. Thanks, Wyatt

My code is as follows:

Main file:

eps = 1.0e-6;

KernelA1(xv)
  for kl = 1:100
      dxv = inv(JacobianA1(xv))*KernelA1(xv);
      xv = xv - dxv/5;
      resA1 = norm(dxv)
      if norm(dxv) < eps
          break;
      end
  end

Jacobian function:

function Jv = JacobianA1(xv)

n = length(xv);
Jv = zeros(n,n);

for ki = 1:n
  fi = KernelA1(xv);

    for kj = 1:n
        if abs(xv(kj)) > 1e-12
            dh = xv(kj)/1e4;
        else
            dh = 1e-6;
        end

        pert = xv(kj) + dh;
        xvnew = [xv(1:kj-1);pert;xv(kj+1:n)];
        fij = KernelA1(xvnew);

        Jv(ki,kj) = (fij(ki)-fi(ki))/dh;
    end
end