%Muller's method to find the roots of the polynomial
%f(x)=x^3-x^2+3x-2
%uses function calls
%ME2450 Eric Pardyjak

clc
clear
itermax=100 %set the maximum number of iterations
x=-1:.1:4;  

%initialize Guesses based on graph to try to obtain root near x=2
%sometimes an imaginary guess works
x0=0
x2=2i
estop = 0.01;
zerr = 100;
x1=(x0+x2)/2    
iter=0;

%Iterate 3 times to try and find the root near x=2
while zerr > estop & iter < itermax

   iter=iter+1                     %iteration number
    
%caluclate the function value at x0,x1,x2 by calling f.m
    
    h0=x1-x0;
    h1=x2-x1;
    delta0=(f(x1) - f(x0))/(x1 - x0);
    delta1=(f(x2) - f(x1))/(x2 - x1);
    
%calculate the polynomial coefficients a,b and c
    a=(delta1 - delta0)/(h1 + h0);
    b=a*h1+delta1;
    c=f(x2);

 %calculate the discriminant
    s = sqrt(b.^2-4*a*c);

    if b<0
        x3=x2-2*c/(b-s);
    else
        x3=x2-2*c/(b+s);
    end

%calculate the parabolic function
f_parab = a*(x-x2).^2 + b*(x-x2) + c;
%Plot the function out graphically
%figure
%plot(x,f(x),x,f_parab,'r',x,0,x0,f(x0),'*',x1,f(x1),'*',x2,f(x2),'*')    


%Method to check error with complex roots
    z3=sqrt(real(x3)^2+imag(x3));
    z2=sqrt(real(x2)^2+imag(x2));
    zerr = (abs(z3-z2)/z3)*100

% update guesses
    x0=x1;
    x1=x2;
    x2=x3
    
end
%right out final root to screen
xroot=x3



% plot(x,f_parab,'r')
% plot(x,0)
% plot(x0,f(x0),'*')
% plot(x1,f(x1),'*')
% plot(x2,f(x2),'*')
