Golden Section Vs Fminbnd Function

clc;
clear;
close all;

L=1000; %feet

Q=20; %gmp


f=@(D) 0.45*L+(0.245*L*(D.^1.5))+325*((44*10^-8*(L*Q.^3/(D.^5))+1.92*(10^-9)*(L*Q^2.68)/(D^4.68)).^(1/2))+61.6*(44*10^-8*(L*Q.^3/(D.^5))+1.92*(10^-9)*(L*Q^2.68)/(D^4.68)).^0925+102

a=.25;
b=6;
error=0.0001; %Stoping Criteria for Golden Section


%=====================

%Start Golden Section

r=(sqrt(5)-1)/2;
x0=a;
x3=b;
x1=a+(1-r)*(b-a);
x2=a+r*(b-a);
y1=feval(f,x1);
y2=feval(f,x2);
i=0
itteration_goldensection=[];
error_goldensection=[];
while (x3-x0)*r>error
    if y1<=y2
        x3=x2;
        x2=x1;
        x1=x0+(1-r)*(x3-x0);
    else
        x0=x1;
        x1=x2;
        x2=x0+r*(x3-x0);
    end
    y1=feval(f,x1);
    y2=feval(f,x2);
    i=i+1;
    itteration_goldensection=[itteration_goldensection i];
    error_goldensection=[error_goldensection abs(feval(f,(x0+x3)/2)-0.425516)];
end
plot(itteration_goldensection,error_goldensection,'r-.');
y_array=[feval(f,x0),y1,y2,feval(f,x3)];
[ymin,index]=min(y_array);
x_array=[x0,x1,x2,x3];
xmin=x_array(index);
disp(sprintf('Golden Search: minimum of f(x) is %10.6f that occurs at x=%10.6f\n',ymin,xmin))
%====================

%MATLAB FUNCTION

[x, feval]=fminbnd(f,.25,6)