function [xmin, fmin, count] = mullersearch(f, df,searchpos,tol)

%figure(3);
count = 0;
x0 = searchpos(1);
x2 = searchpos(2);
x1 = (x2-x0)/2;
x3 = 0;

ea = realmax;
left = x0;
right = x2;
verbose = true;

a = 0;
b = 0;
c = 0;
figure(3);
set(gcf,'Position',[7 7 700 1000]);
hold on
[miny, maxy] = drawplot(f,df, a,b,c, x0, x1, x2, count);
placelabel(left,0,count);
placelabel(right,0,count);

export_fig(gcf, '1.jpg', '-transparent', '-r300');

%iterative computation

while(ea > tol)

    count = count+1;
    %1. compute for h0,h1,d0,d1

    fx0 = df(x0);
    fx1 = df(x1);
    fx2 = df(x2);

    h0 = x1-x0;
    h1 = x2-x1;
    d0 = ( fx1-fx0 ) / h0;
    d1 = ( fx2-fx1 ) / h1;

    %2. compute for a,b,c
    a = (d1-d0) / (h1-h0);
    b = a*h1 + d1;
    c = fx2;
    %if count == 3
    %end

    %3. check for sign
    temp = real(sqrt( b*b - 4*a*c ));
    checkpos = b + temp;%holder for |b+sqrt(b^2-4ac)|
    checkneg = b - temp;%holder for |b-sqrt(b^2-4ac)|

    %4. compute for x3, using form according to sign
    if abs(checkneg) < abs(checkpos)
        x3 = x2 - (2*c) / checkpos;
    else
        x3 = x2 - (2*c) / checkneg;
    end

    drawplot(f, df,  a, b, c, x0, x1, x2, count);
    if count <= 10
        export_fig(gcf,[num2str(count),'.jpg'],'-transparent','-r300');
    end

    %5. compute approximate error
    ea = abs((x3-x2)/x3);

    %6. if verbose, print values
    if verbose == true

        % print: count, x0, x1, x2, x3, ea
        sprintf('Iteration: %i' , count)
        sprintf('x0 = %f', x0)
        sprintf('x1 = %f', x1)
        sprintf('x2 = %f', x2)
        sprintf('Root Estimate: x3 = %f', x3)
        sprintf('Approximate Relative Error: %f', ea)
    end


    %7. adjust values for next iteration
    x0 = x1;
    x1 = x2;
    x2 = x3;

end
%worst case scenario: count maxed out
sprintf('Maximum iterations exhausted')
xmin = x3;
fmin = df(x3);
end

function [miny, maxy] = drawplot(f, df, A, B, C, x1, x2, x3, count)
h = 0.05;
asixLeft = - 20;
asixRight = 20;
left = asixLeft -100;
right = asixLeft +100;

x = left:h:right;
y = f(x);
miny = min(y) - 10;
maxy = max(y) + 10;
dy = df(x);
mindy = min(dy) - 10;
maxdy = max(dy) + 10;

colp = hsv2rgb([rand(), 1, 0.5+0.5*rand()]);
hyper = @(x) A*x.^2 + B*x + C;
yh = hyper(x);
col = hsv2rgb([rand(), 1, 0.5+0.5*rand()]);

newplot;
hold on
xlabel('\itx');
ylabel('\ity');
scatter([left right],[f(left), f(right)],'Marker','o','MarkerFaceColor',colp,'MarkerEdgeColor',colp);

newplot;
hold on
col = hsv2rgb([rand(), 1, 0.5+0.5*rand()]);
subplot(2,1,1);
line([left right],[0 0],'Color','k','LineWidth',1); %axis x
plot(x,yh,'Color',col, 'LineWidth',0.5);
y1 = hyper(x1);

line([x1 x1],[0 y1],'Marker','s','Color',col,'LineWidth',1,'MarkerSize',4);
y2 = hyper(x2);
line([x2 x2],[0 y2],'Marker','s','Color',col,'LineWidth',1,'MarkerSize',4);
y3 = hyper(x3);
line([x3 x3],[0 y3],'Marker','s','Color',col,'LineWidth',1,'MarkerSize',4);
plot(x,dy,'Color',colp, 'LineWidth',2);
placelabel(x1,y1,count);
placelabel(x2,y2,count);
placelabel(x3,y3,count);
axis([asixLeft asixRight -20000 20000]);

newplot;
hold on
subplot(2,1,2);
line([left right],[0 0],'Color','k','LineWidth',1); %axis x
y1 = f(x1);
line([x1 x1],[0 y1],'Marker','s','Color',col,'LineWidth',1,'MarkerSize',4);
y2 = f(x2);
line([x2 x2],[0 y2],'Marker','s','Color',col,'LineWidth',1,'MarkerSize',4);
y3 = f(x3);
line([x3 x3],[0 y3],'Marker','s','Color',col,'LineWidth',1,'MarkerSize',4);
plot(x,y,'Color',colp, 'LineWidth',2);
axis([asixLeft asixRight -20000 20000]);
end


function placelabel(x,y,iternumber)
if iternumber <=10
    text(x, 0, num2str(iternumber), 'BackgroundColor', 'white');
end
end