%Program 1. Bisection Method
%Computes approximate solution of f(x)=0
%Input: inline function f; a,b such that f(a)*f(b)<0, 
%     and iterartion number n positive integer
%Output: Bracketing interval [an,bn] and approximate  solution cn.
function [an,cn,bn] = bisect(f,a,b,n)
if sign(f(a))*sign(f(b)) >= 0
  error('f(a)f(b)<0 not satisfied!') %ceases execution
end
an=a; bn=b; c=(a+b)/2;
fa=f(a);
fb=f(b);
k = 1;
while k<n
  c=(a+b)/2;
  fc=f(c);
  if fc == 0              %c is a solution, done
    break
  end
  if sign(fc)*sign(fa)<0  %a and c make the new interval
    b=c;fb=fc; 
  else                    %c and b make the new interval
    a=c;fa=fc;
  end
  k=k+1;
end
an=a;
bn=b;
cn=(an+bn)/2;

___________________________________________________________________________
%Command Window:
%Example 1. Solving f(x)=x^3+4x^2-10=0
clear;
format long
f=inline('x^3+4*x^2-10','x');
A=zeros(20:4);
for i=1:20;
[x,y,u,w]=bisect1(f,1,2,i);
A(i,1)=x; %left end point
A(i,2)=y; %mid point
A(i,3)=u; %right end point
A(i,4)=w; %function value
 end
disp(' ')    %  a blank line
disp(' n         an                  cn                  bn          f(cn) ')
disp('-----------------------------------------------------------------------')
for k=1:20
  disp(sprintf('%2.0f    %7.10e    %7.10e   %7.10e  %7.10e',k,A(k,1),A(k,2),A(k,3),A(k,4)))
end
 
 n         an                  cn                  bn          f(cn) 
-----------------------------------------------------------------------
 1    0.0000000000e+000    1.0000000000e+000   2.0000000000e+000  -5.0000000000e+000
 2    1.0000000000e+000    1.5000000000e+000   2.0000000000e+000  2.3750000000e+000
 3    1.0000000000e+000    1.2500000000e+000   1.5000000000e+000  -1.7968750000e+000
 4    1.2500000000e+000    1.3750000000e+000   1.5000000000e+000  1.6210937500e-001
 5    1.2500000000e+000    1.3125000000e+000   1.3750000000e+000  -8.4838867188e-001
 6    1.3125000000e+000    1.3437500000e+000   1.3750000000e+000  -3.5098266602e-001
 7    1.3437500000e+000    1.3593750000e+000   1.3750000000e+000  -9.6408843994e-002
 8    1.3593750000e+000    1.3671875000e+000   1.3750000000e+000  3.2355785370e-002
 9    1.3593750000e+000    1.3632812500e+000   1.3671875000e+000  -3.2149970531e-002
10    1.3632812500e+000    1.3652343750e+000   1.3671875000e+000  7.2024762630e-005
11    1.3632812500e+000    1.3642578125e+000   1.3652343750e+000  -1.6046690755e-002
12    1.3642578125e+000    1.3647460938e+000   1.3652343750e+000  -7.9892628128e-003
13    1.3647460938e+000    1.3649902344e+000   1.3652343750e+000  -3.9591015229e-003
14    1.3649902344e+000    1.3651123047e+000   1.3652343750e+000  -1.9436590101e-003
15    1.3651123047e+000    1.3651733398e+000   1.3652343750e+000  -9.3584728188e-004
16    1.3651733398e+000    1.3652038574e+000   1.3652343750e+000  -4.3191879925e-004
17    1.3652038574e+000    1.3652191162e+000   1.3652343750e+000  -1.7994890323e-004
18    1.3652191162e+000    1.3652267456e+000   1.3652343750e+000  -5.3962541529e-005
19    1.3652267456e+000    1.3652305603e+000   1.3652343750e+000  9.0309927430e-006
20    1.3652267456e+000    1.3652286530e+000   1.3652305603e+000  -2.2465803845e-005
% Example 2. Solving f(x)=xsin(x)-1=0.
% Finding a nad b, by plotting

 format long
 f=inline('x*sin(x)-1','x');
A=zeros(20:4);
 for i=1:20;
[x,y,u,w]=bisect1(f,0,2,i);
A(i,1)=x; %left end point
A(i,2)=y; %mid point
A(i,3)=u; %right end point
A(i,4)=w; %function value
 end
disp(' ')    %  a blank line
disp(' n         an                  cn                  bn          f(cn) ')
disp('-----------------------------------------------------------------------')
for k=1:20
  disp(sprintf('%2.0f    %7.10e    %7.10e   %7.10e  %7.10e',k,A(k,1),A(k,2),A(k,3),A(k,4)))
end
 
 n         an                  cn                  bn               f(cn) 
------------------------------------------------------------------------------
 1    0.0000000000e+000    1.0000000000e+000   2.0000000000e+000  -1.5852901519e-001
 2    1.0000000000e+000    1.5000000000e+000   2.0000000000e+000  4.9624247991e-001
 3    1.0000000000e+000    1.2500000000e+000   1.5000000000e+000  1.8623077419e-001
 4    1.0000000000e+000    1.1250000000e+000   1.2500000000e+000  1.5051043361e-002
 5    1.0000000000e+000    1.0625000000e+000   1.1250000000e+000  -7.1826631385e-002
 6    1.0625000000e+000    1.0937500000e+000   1.1250000000e+000  -2.8361722663e-002
 7    1.0937500000e+000    1.1093750000e+000   1.1250000000e+000  -6.6427748603e-003
 8    1.1093750000e+000    1.1171875000e+000   1.1250000000e+000  4.2080340106e-003
 9    1.1093750000e+000    1.1132812500e+000   1.1171875000e+000  -1.2164904180e-003
10    1.1132812500e+000    1.1152343750e+000   1.1171875000e+000  1.4960036578e-003
11    1.1132812500e+000    1.1142578125e+000   1.1152343750e+000  1.3981310236e-004
12    1.1132812500e+000    1.1137695313e+000   1.1142578125e+000  -5.3832472254e-004
13    1.1137695313e+000    1.1140136719e+000   1.1142578125e+000  -1.9925230310e-004
14    1.1140136719e+000    1.1141357422e+000   1.1142578125e+000  -2.9718720730e-005
15    1.1141357422e+000    1.1141967773e+000   1.1142578125e+000  5.5047411086e-005
16    1.1141357422e+000    1.1141662598e+000   1.1141967773e+000  1.2664400201e-005
17    1.1141357422e+000    1.1141510010e+000   1.1141662598e+000  -8.5271465147e-006
18    1.1141510010e+000    1.1141586304e+000   1.1141662598e+000  2.0686302813e-006
19    1.1141510010e+000    1.1141548157e+000   1.1141586304e+000  -3.2292572574e-006
20    1.1141548157e+000    1.1141567230e+000   1.1141586304e+000  -5.8031327321e-007
% Multiple root: Solve f(x)=(x-1/3)^3=x^3-x^2+x/3-1/27
clear;
format long
f=inline('x^3-x^2+x/3-1/27','x');
A=zeros(25:4);
for j=1:25;
[x,y,u,w]=bisect1(f,0,1,j);
A(j,1)=x; %left end point
A(j,2)=y; %mid point
A(j,3)=u; %right end point
A(j,4)=w; %function value
 end
disp(' ')    %  a blank line
disp(' n         an                  cn                  bn          f(cn) ')
disp('-----------------------------------------------------------------------')
for k=1:25
  disp(sprintf('%2.0f    %7.10e    %7.10e   %7.10e  %7.10e',k,A(k,1),A(k,2),A(k,3),A(k,4)))
end
 
 n         an                  cn                  bn          f(cn) 
-----------------------------------------------------------------------
 1    0.0000000000e+000    5.0000000000e-001   1.0000000000e+000  4.6296296296e-003
 2    0.0000000000e+000    2.5000000000e-001   5.0000000000e-001  -5.7870370370e-004
 3    2.5000000000e-001    3.7500000000e-001   5.0000000000e-001  7.2337962963e-005
 4    2.5000000000e-001    3.1250000000e-001   3.7500000000e-001  -9.0422453704e-006
 5    3.1250000000e-001    3.4375000000e-001   3.7500000000e-001  1.1302806713e-006
 6    3.1250000000e-001    3.2812500000e-001   3.4375000000e-001  -1.4128508391e-007
 7    3.2812500000e-001    3.3593750000e-001   3.4375000000e-001  1.7660635496e-008
 8    3.2812500000e-001    3.3203125000e-001   3.3593750000e-001  -2.2075794387e-009
 9    3.3203125000e-001    3.3398437500e-001   3.3593750000e-001  2.7594743157e-010
10    3.3203125000e-001    3.3300781250e-001   3.3398437500e-001  -3.4493422008e-011
11    3.3300781250e-001    3.3349609375e-001   3.3398437500e-001  4.3116760162e-012
12    3.3300781250e-001    3.3325195313e-001   3.3349609375e-001  -5.3895776730e-013
13    3.3325195313e-001    3.3337402344e-001   3.3349609375e-001  6.7376659807e-014
14    3.3325195313e-001    3.3331298828e-001   3.3337402344e-001  -8.4238171993e-015
15    3.3331298828e-001    3.3334350586e-001   3.3337402344e-001  1.0547118734e-015
16    3.3331298828e-001    3.3332824707e-001   3.3334350586e-001  -1.2490009027e-016
17    3.3332824707e-001    3.3333587646e-001   3.3334350586e-001  1.3877787808e-017
18    3.3332824707e-001    3.3333206177e-001   3.3333587646e-001  0.0000000000e+000
19    3.3332824707e-001    3.3333206177e-001   3.3333587646e-001  0.0000000000e+000
20    3.3332824707e-001    3.3333206177e-001   3.3333587646e-001  0.0000000000e+000
21    3.3332824707e-001    3.3333206177e-001   3.3333587646e-001  0.0000000000e+000
22    3.3332824707e-001    3.3333206177e-001   3.3333587646e-001  0.0000000000e+000
23    3.3332824707e-001    3.3333206177e-001   3.3333587646e-001  0.0000000000e+000
24    3.3332824707e-001    3.3333206177e-001   3.3333587646e-001  0.0000000000e+000
25    3.3332824707e-001    3.3333206177e-001   3.3333587646e-001  0.0000000000e+000