% This MATLAB code implements the Newton method to solve for roots of the
% input equation myfun with derivative myfunprime.
% Call as bisection(a,b,'myfun',tol)
function[approx]=bisection(a,b,myfun,tol)
% Set up iteration control
maxiter=10000; i = 1;
% Set up initial approx
approx = (a+b)/2;
while abs(feval(myfun,approx))>=tol && i<=maxiter;
    c=(a+b)/2;
    % Implement the Bisection method
    %determine where point c lies and update interval bounds accordingly
    if feval(myfun,c)*feval(myfun,a)>0; 
        a=c;
    else;
        b=c;
    end;
    approx=(a+b)/2;
    % Update iteration count
    i=i+1;
    % Print a message if maximum iteraton is reached:
    if i==maxiter;
        fprintf('Max number of %d iterations reached; root may not\n',maxiter);
        fprintf('have been found to within tolerance\n');
    end;
end;
fprintf('%d iterations',i);