clear       // clear all previous variables



function bisection(a,b,tolerance,maxstep)
    // The function takes four inputs: the start and end points of the original interval, 
    // 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 on the requested interval, to see where the root might be
    clf
    x = [ a: 0.1 :b ]
    plot(x,f(x))
    g = gca();  
    g.x_location = "origin"; 
    g.y_location = "origin"; 
    
    // Now for the method itself:
    step = 1        // Keep track of which step we're on
    leftY = f(a)    // Remember first left endpoint before we change a

    converge = %F   // Assume the process won't work, change later if it does
     
    while (step <= maxstep)     // Keep repeating, unless maximum steps reached   
        p = a + (b-a)/2         // midpoint of interval
 
         
        if(f(p) == 0 | abs(f(p)) < tolerance) then
            disp('Step '+string(step))
            disp('     Interval: ['+string(a)+', '+string(b)+']')
            disp('     Midpoint, p: '+string(p))
            disp('     f(p): '+string(f(p)))

            converge = %T       // we found the root, or close enough, so set converge to True
            break
        else
            // no root found yet, print step info
            disp('Step '+string(step))
            disp('     Interval: ['+string(a)+', '+string(b)+']')
            disp('     Midpoint, p: '+string(p))
            disp('     f(p): '+string(f(p)))
            // add a marker to the axis
            plot(p,f(p),'d')
        end
         
        // increase step counter  
        step = step + 1
         
        // see which half of interval the function changes sign, and use as next interval
        if(leftY*f(p) > 0) then
            a = p
        else
            b = p
        end
    end

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

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

bisection( -5, 5, 0.0001, 100 )
