parabolic.m

function [E4,x4] = parabolic(x_l,x_u)
epsilon = .039;
sigma = 2.934;
fun = @(x) 4*epsilon*((sigma./x).^12-(sigma./x).^6);
x = 3.2933;
x_l = 3.22;
x_u = 3.45;

x1 = x_l;
x2 = mean([x_l,x_u]);
x3 = x_u;

f1 = fun(x1);
f2 = fun(x2);
f3 = fun(x3);
p = polyfit([x1,x2,x3],[f1,f2,f3],2);
x_fit = linspace(x1,x3,20);
y_fit = polyval(p,x_fit);

plot(x,fun(x),x_fit,y_fit,[x1,x2,x3],[f1,f2,f3],'o')
hold on
if f2<f1 && f2<f3
    x4=x2-0.5*((x2-x1)^2*(f2-f3)-(x2-x3)^2*(f2-f1))/((x2-x1)*(f2-f3)-(x2-x3)*(f2-f1));
    f4=fun(x4);

    if x4>x2
        plot(x4,f4,'*',[x1,x2],[f1,f2])
        x1=x2;
        f1=f2;
    else
        plot(x4,f4,'*',[x3,x2],[f3,f2])
        x3=x2;
        f3=f2;
    end
    x2=x4; f2=f4;
else
    %error('no minimum in bracket')
    f4 = 'no';
end
hold off
E4=f4