clear       // clear all previous variables

function secant(x0,x1,tolerance,maxstep)
    // The function takes four inputs: the two starting guesses of the root, 
    // the required tolerance, and the max number of steps at which to terminate in case
    // the method does not converge before then (always a good idea)
    
    // Draw a picture of the function near the starting guesses:
    clf
    x = [ x0-2 : 0.1 : x1+2 ]
    plot(x,f(x))
    g = gca();  
    g.x_location = "origin"; 
    g.y_location = "origin"; 
    plot(x0,f(x0),'s')
    plot(x1,f(x1),'s')
    
    // Now for the method itself:
    step = 1        // Keep track of which step we're on
 
    converge = %F   // Assume the process won't work, change later if it does
     
    while (step <= maxstep)     // Keep repeating, unless maximum steps reached   
        p = x1 - f(x1)*((x1-x0)/(f(x1)-f(x0)))         // get next approximation
        
         
        if(f(p) == 0 | abs(f(p)) < tolerance) then
            disp('Step '+string(step))
            disp('     x_i: '+string(x0) )
            disp('     x_i+1: '+string(x1) )
            disp('     Next guess, x_i+2: '+string(p))
            disp('     f(x_i+2): '+string(f(p)))

            converge = %T
            break
        else
            // no root found yet, print step info
            disp('Step '+string(step))
            disp('     x_i: '+string(x0) )
            disp('     x_i+1: '+string(x1) )
            disp('     Next guess, x_i+2: '+string(p))
            disp('     f(x_i+2): '+string(f(p)))
            // add a marker to the axis
            
            plot(p,f(p),'d')
        end
         
        // increase step counter  
        step = step + 1
         
        // redefine starting pair:
        x0 = x1;
        x1 = p;
        
    end

    // print summary
    if (converge) then
        disp('Approximate root: '+string(p)+' found in '+string(step)+' steps')
    else
        disp('Root not found! Try different starting points or adjust tolerance')
    end
     
endfunction

function f = f(x)
    f = x^2 - 4        // change to your favorite function
endfunction

secant( 1.5, 2.8, 0.0001, 100 )
