The following code is my attempt of employing the Broydens method on the function $f(x,y)=(e^{xy}-y^2-2,\cos(x+y)+\frac{1}{2})$. Where the first Matrix is the Jacobian, then it gets updated using an estimate. So $k[l,k]$ computes the first 12 iterations. However, this program is extremely slow and I am not sure why. Could anyone point out why it is so slow and how I could improve it. If you are not familiar https://en.wikipedia.org/wiki/Broyden%27s_method described the Broyden Method nicely.

EDIT: I found some optimizaiton to my code. I removed the redundant $N$ as per a suggestion in the comments. The optiumization I found is to deal with numerical matricies form the get go. The original code did everything for arbitrary $x,y$

The following code is my attempt of employing Broyden's method for root finding on the function $f(x,y)=(e^{xy}-y^2-2,\cos(x+y)+\frac{1}{2})$. Where the first matrix is the Jacobian, then it gets updated using an estimate. So $k[l,k]$ computes the first 12 iterations. However, this program is extremely slow and I am not sure why. Could anyone point out why it is so slow and how I could improve it?

EDIT: I found some optimization to my code. I removed the redundant $N$ as per a suggestion in the comments. The optimization I found is to deal with numerical matrices form the get go. The original code did everything for arbitrary $x,y$

I would appriciate any input.

I would appreciate any input.

EDIT: I found some optimizaiton to my code. I removed the redundant $N$ as per a suggestion in the comments. The optiumization I found is to deal with numerical matricies form the get go. The original code did everything for arbitrary $x,y$

k[l_, k_] := Module[{m}, 
  f[{x_, y_}] := {Exp[x y] - y^2 - 2, Cos[x + y] + 1/2};
  h = {Exp[x y] - y^2 - 2, Cos[x + y] + 1/2};
  b = {x, y};
  
  m[0] := N[{{l}, {k}}];
  
  J[0] := 
   N[D[h, {b}]] /. {x -> Part[Part[m[0], 1], 1], 
     y -> Part[Part[m[0], 2], 1]};
  m[n_] := m[n] = m[n - 1] - LinearSolve[J[n - 1], f[m[n - 1]]];
  J[n_] := 
   J[n] = J[
      n - 1] + ((f[m[n]] - f[m[n - 1]] - J[n - 1].(m[n] - m[n - 1]))/
        Norm[m[n] - m[n - 1]]^2).Transpose[m[n] - m[n - 1]];
  Table[m[a], {a, 1, 15}]]
```