channelflow.org

A bare-bones Newton search code for solving the nonlinear equation f(x) = 0.

function x = newtonsearch(f,x)
% return x that solves f(x) == 0, using input argument x as initial guess

eps = 1e-14;  % tolerance 1e-14 == 1 x 10^(-14)
dx  = 1e-07;  % increment for approximating slope df/dx

fx = f(x);

while (abs(fx) > eps) 
 
  % approximate slope f'(x) 
  fx_dx = f(x+dx);
  fprime = (fx_dx - fx)/dx; 
  
  % solve for newton step delta_x
  delta_x = -fx/fprime;
  
  % set new guess for x to x + delta_x and compute new f(x)
  x = x + delta_x;
  fx = f(x);
  
end