function [r]=Initial(x_guess,p) A=GetJacobianInv2(x_guess,p); [r] = Fcn2(x_guess,p); err = norm(r.out); err1 = err; inner_ct = 0; x_new = x_guess; while (err > p.tol) && (inner_ct < p.max_inner_ct) inner_ct = inner_ct + 1; smaller = 0; j = 0; f0 = r.out; d = -(A*r.out); d2(p.n+1,1)=0; d2(1:p.n,1)=d(1:p.n,1); for j = 0:p.max_steps-1 x_new = x_guess + (2^-j)*d2; [r]=Fcn2(x_new,p); err1 = norm(r.out); if (err1 < err) break; end end w = r.out - f0; q = x_new - x_guess; q2(1:p.n,1)=q(1:p.n,1); qA = (q2'*A); denom = qA*w; A_update = (A*w - q2)*qA/denom; %%%%%%%%%%Sherman-Morrison Formula%%%%%%%%%%% A = A - A_update; if mod(inner_ct, 20) == 0 [1 err1] A = GetJacobianInv2(x_new, p); inner_ct = inner_ct + 5; end if abs(denom) < 1e-20 || norm(q) < 1e-15 || ... norm(A_update) > 1e8 [2 denom norm(q) norm(A_update) err1] A = GetJacobianInv2(x_new, p); inner_ct = inner_ct + 20; end err = err1; x_guess = x_new; end