Octave code for Muller's method

% Muller's method
function rtn = mullers(fx,x2,x0,x1,n)

X =[];

for i = 1:n
  x = x2; f2 = eval(fx);
  x = x0; f0 = eval(fx);
  x = x1; f1 = eval(fx);

  h1 = x1 - x0;
  h2 = x0 - x2;
  gamma = h2 / h1;

  c = f0;
  a = (gamma*f1 - f0*(1+gamma) + f2) / (gamma * h1^2 * (1+gamma));
  b = (f1 - f0 - a*h1^2) / h1;

  if b>0
    root = x0 - 2*c / (b + sqrt(b^2 - 4*a*c));
  elseif
    root = x0 - 2*c / (b - sqrt(b^2 - 4*a*c));
  endif

  if root > x0
    x2 = x0;
    if root > x1
      x0 = x1; x1 = root;
    else
      x0 = root;
    endif
  elseif
    x1 = x0;
    if root < x2
      x0 = x2; x2 = root;
    else
      x0 = root;
    endif
  endif

  X = [X; i, root];
endfor

disp (X)

endfunction

octave:1> fx = '3*x+sin(x)-exp(x)'
fx = 3*x+sin(x)-exp(x)
octave:2> mullers(fx,0.0,0.5,1.0,3)
warning: suggest parenthesis around assignment used as truth value
warning: suggest parenthesis around assignment used as truth value
   1.00000   0.35491
   2.00000   0.36046
   3.00000   0.36042
octave:3>