function [err,n5ntvalue,n1,n2,n3,n4,n5,n6,Itterations] = BisectionErrSolve(a,b,...
    max_iterations,tolerance,pressure,temperature,fuelequiv,carbon,...
    oxygen,hydrogen,molar_air_stoich,gibbseng)

% This function uses a while loop for the bisection method to root solve 
% for n5nt using the 'nEquations' function to calculate n5nt and error. 

[kco2,kh2o] = GibbsEngCalcs(temperature,gibbseng);

iteration = 0;

while abs((b - a) / 2) > tolerance && iteration < max_iterations
    iteration = iteration + 1;
    c = (a + b) / 2;
    
    % Evaluate error at midpoint
    f_c = nEquations(c, pressure,fuelequiv,...
    carbon,oxygen,hydrogen,molar_air_stoich,kco2,kh2o);
    
    f_a = nEquations(a, pressure,fuelequiv,...
    carbon,oxygen,hydrogen,molar_air_stoich,kco2,kh2o);
    
    if f_c == 0
        % Exact root found
        break;
    elseif sign(f_c) * sign(f_a) < 0
        b = c; % Root is in [a, c]
    else
        a = c; % Root is in [c, b]
    end
    Itterations = iteration;
end

n5ntvalue = (a + b) / 2; % Approximate root

[err,n5ntvalue,n1,n2,n3,n4,n5,n6] = nEquations(n5ntvalue, pressure,fuelequiv,...
    carbon,oxygen,hydrogen,molar_air_stoich,kco2,kh2o);

% fprintf(['Bisection method converged n5ntvalue: %.6f after %d' ...
%     ' iterations at error Size: %.14f\n'], n5ntvalue, iteration,err);
% 
% fprintf(['Values Calculated: n5nt = %.3f, n1 = %.3f, n2 = %.3f,'...
% ' n3 = %.3f, n4 = %.3f, n5 = %.3f, n6 = %.3f\n\n'],n5ntvalue,n1,n2,n3,n4,n5,n6);
end