clear, clc

format short g
% solution to exercise 4.2

% (1) to find three intervals with length less than 0.1 containing all
% three solutions

% use the following plot to guess the three solution locations

x=-10:0.1:10;
y=feval(@exercise42,x); 

plot(x,y); grid; 

disp('It can be oberved that the three solutions roughly locate in interval [-0.2,0], [0,0.2] and [8.9,9.1]')

% define the intervals 

a1=-0.4; b1=0; 
a2=0; b2=0.4; 
a3=8.8; b3=9.2; 
k=0;

while (b1-a1>0.1) | (b2-a2>0.1) | (b3-a3>0.1)
    m1=(a1+b1)/2;
    m2=(a2+b2)/2;
    m3=(a3+b3)/2;
    fa1=feval(@exercise42,a1); 
    fb1=feval(@exercise42,b1); 
    fm1=feval(@exercise42,m1);
    if fa1*fm1 <0
        b1=m1;
    else
        a1=m1;
    end
    
    fa2=feval(@exercise42,a2); 
    fb2=feval(@exercise42,b2); 
    fm2=feval(@exercise42,m2);
    if fa2*fm2 <0
        b2=m2;
    else
        a2=m2;
    end
    
    fa3=feval(@exercise42,a3); 
    fb3=feval(@exercise42,b3); 
    fm3=feval(@exercise42,m3);
    if fa3*fm3 <0
        b3=m3;
    else
        a3=m3;
    end
    k=k+1;
end

disp('the intervals are')
a1
b1
a2
b2
a3
b3

disp('the no.of iteration is'); k


%(2) verification is skip...

%(3) determine which expression will converge to which solution

% first try starts at a1 point
p11old=a1;
p12old=a1;
p13old=a1;
x=zeros(10,3);

for k=1:10
    p11=feval(@exfunc3,p11old);
    p12=feval(@exfunc4,p12old);
    p13=feval(@exfunc5,p13old);
    x(k,:)=[p11 p12 p13];
    p11old=p11; 
    p12old=p12;
    p13old=p13;
end

disp('  #Iterations    eq1       eq2   eq3')
disp([[1:10]' x])

disp('starting from a1=-0.1, eq1 converges to 0.10541, eq3 converges to -0.095345, while eq2 diverges')

disp('press any key for continuing...')
pause

% second try starts at a2 point
p11old=a2;
p12old=a2;
p13old=a2;
x=zeros(10,3);

for k=1:10
    p11=feval(@exfunc3,p11old);
    p12=feval(@exfunc4,p12old);
    p13=feval(@exfunc5,p13old);
    x(k,:)=[p11 p12 p13];
    p11old=p11; 
    p12old=p12;
    p13old=p13;
end

disp('  #Iterations    eq1       eq2   eq3')
disp([[1:10]' x])

disp('starting from a2=0.1, eq1 converges to 0.10541, eq3 converges to -0.095345, while eq2 diverges')

disp('press any key for continuing...')
pause

% third try starts at a2 point
p11old=a3;
p12old=a3;
p13old=a3;
x=zeros(10,3);

for k=1:10
    p11=feval(@exfunc3,p11old);
    p12=feval(@exfunc4,p12old);
    p13=feval(@exfunc5,p13old);
    x(k,:)=[p11 p12 p13];
    p11old=p11; 
    p12old=p12;
    p13old=p13;
end

disp('  #Iterations    eq1       eq2   eq3')
disp([[1:10]' x])

disp('starting from a3=8.9, eq1 converges to 0.10541, eq3 converges to -0.095345, while eq2 converges to 8.9995 ')

disp('press any key for continuing...')
pause


%(4) solution up to tolerance 10e-6

format long g
p11old=a3;
p12old=a3;
p13old=a3;
p11=feval(@exfunc3,p11old);
p12=feval(@exfunc4,p12old);
p13=feval(@exfunc5,p13old);

while abs(p11-p11old)>10e-6 | abs(p12-p12old)>10e-6 | abs(p13-p13old)>10e-6
    p11old=p11; 
    p12old=p12;
    p13old=p13;
    p11=feval(@exfunc3,p11old);
    p12=feval(@exfunc4,p12old);
    p13=feval(@exfunc5,p13old);  
end

disp('  the final solutions with tolerance 10e-6 are')
p11
p12
p13

disp('The exercise is over. ')

