function [x_sol, f_at_x_sol, N_iterations] = bisection(f, xn, xp, eps_f, eps_x)
% solving f(x)=0 with bisection method
% 
%   f is the function handle to the desired function,
%   xn and xp borders of search,
%   f(xn)<0 and f(xp)>0 required,
%   eps_f defines how close f(x) should be to zero,
%   eps_x defines uncertainty of solution x

%  We need to sanitize the inputs but I skip it for simplicity

	xg=(xp+xn)/2; % initial guess
	fg=f(xg);     % initial function evaluation

	N_iterations=0;      % initialization of the counter

	while ( (abs(fg) > eps_f) & (abs(xg-xp)>eps_x) )
		if (fg>0) 
			xp=xg;
		else
			xn=xg;
		end
		xg=(xp+xn)/2; % update guess
		fg=f(xg);     % update function evaluation

		N_iterations=N_iterations+1;
	end
	x_sol=xg; % solution is ready
	f_at_x_sol=fg;
end
