function trabalho01_ex2a()
    molas = 100;
    variaveis = 2 * (molas + 1) - 4;
    x0 = zeros(1, variaveis);
    for i=1:2:variaveis
        if i <= variaveis/2
            percentual = i / (variaveis/2);
            x0(i) = 0.01 * percentual;
            x0(i+1) = -0.1 * percentual;
        else
            percentual = (variaveis-i) / (variaveis/2);
            x0(i) = 0.01 * percentual;
            x0(i+1) = -0.1 * percentual;
        end
    end

    for i=1:2:variaveis
        u(floor(i/2)+1) = x0(i);
        v(floor(i/2)+1) = x0(i+1);
    end

    x0 = [u,v];

    %x0 = [0.01, -0.1];
    comprimento1 = 30;
    comprimento2 = 30;
    setComprimento(comprimento1, comprimento2, molas);
    setRho(8, 16);
    setEA(27000, 18000);
    %setRho(8, 8);
    %setEA(27000, 27000);

    global g_l;
    global g_k;
    for k = 1:molas
        g_l(k) = getL(k, molas);
        g_k(k) = getK(k, molas);
    end

    global g_w;
    nos_internos = molas - 1;
    for k = 1:nos_internos
        mola_antes = k;
        mola_depois = k + 1;
        P_antes = getP(mola_antes, molas);
        P_depois = getP(mola_depois, molas);
        g_w(k) = (P_antes/2) + (P_depois/2);
    end

    [x1,iter] = univariante(x0, @f)

    x = zeros(1,molas+1);
    y = zeros(1,molas+1);
    for i=1:molas
        xi = (comprimento1+comprimento2) / molas * i;
        x(i+1) = xi;
        if i < molas
            y(i+1) = -x1(molas + i - 1);
        end
    end
    plot(x, y);
    pause;
end

function setComprimento(comprimento1, comprimento2, num_molas)
    global L1;
    global L2;
    L1 = comprimento1 / (num_molas / 2);
    L2 = comprimento2 / (num_molas / 2);
end

function setRho(densidade1, densidade2)
    global rho1;
    global rho2;
    rho1 = densidade1;
    rho2 = densidade2;
end

function setEA(modulo1, modulo2)
    global EA1;
    global EA2;
    EA1 = modulo1;
    EA2 = modulo2;
end

function y = getL(indice, num_mola)
    global L1;
    global L2;
    if indice <= (num_mola)/2
        y = L1;
    else
        y = L2;
    end
end

function y = getP(indice, num_mola)
    global rho1;
    global rho2;
    global L1;
    global L2;
    if indice <= (num_mola)/2
        y = rho1 * L1;
    else
        y = rho2 * L2;
    end
end

function y = getK(indice,num_mola)
    global EA1;
    global EA2;
    global L1;
    global L2;
    if indice <= (num_mola/2)
        y = EA1/L1;
    else
        y = EA2/L2;
    end
end

function y = g(L)
    global g_d;
    global g_x;
    global g_f;
    y = g_f(g_d*L + g_x);
end

function [u,iter] = univariante(x0,f)
    global g_d;
    global g_x;
    global g_f;
    g_f = f;
    g_x = x0;
    tol = 1e-3;
    n = length(x0);

    g_d = zeros(1,n);
    for k = 1:5000
        xi = g_x;
        for i = 1:n
            g_d(i) = 1;
            [a1,b1] = passoconstante(0, @g, 100, 1);
            c1 = bissecao(a1, b1, 1e-6, tol/10, @g);
            x1 = g_x + g_d*c1;

            g_x = x1;
            g_d(i) = 0;
        end

        norm(x1 - xi)
        if norm(x1 - xi) <= tol
            break;
        end
    end

    u = x1;
    iter = k;
end

function y = f(uv)
    global g_k;
    global g_l;
    tamanho_u = length(uv)/2;

    u = uv(1:tamanho_u);
    v = uv(tamanho_u+1:end);
    num_molas = tamanho_u + 1;

    du = u - [0, u(1:end-1)];
    dv = v - [0, v(1:end-1)];
    du(num_molas) = -u(num_molas-1);
    dv(num_molas) = v(num_molas-1);

    U = 0.5 * sum(g_k .* ((hypot(g_l + du, dv) - g_l).^2));

    global g_w;
    V = sum(g_w .* v);

    y = U - V;
end

function [y,w] = passoconstante(x0,f,da,d)
   y0 = f(x0);
    for k = 1:5000
        x1 = x0 + da*d;
        y1 = f(x1);
        if k == 1
            if y1 > y0
                d = -d;
                continue
            elseif y1 == y0
                y = min([x0 - da*d,x0 + da*d]);
                w = max([x0 - da*d,x0 + da*d]);
                break
            end
        else
            if y1 >= y0
                y = min([x0 - da*d,x0 + da*d]);
                w = max([x0 - da*d,x0 + da*d]);
                break
            end
        end
        x0 = x1;
        y0 = y1;
    end
end

function y = bissecao(aL,aU,eps,tol,f)
    for k = 1:5000
        aM = (aL + aU)/2;
        a1 = aM - eps;
        a2 = aM + eps;
        if abs(aL - aU) <= tol
            y = aM;
            break
        elseif f(a2) > f(a1)
            aU = aM;
        elseif f(a1) > f(a2)
            aL = aM;
        else
            y = aM;
            break
        end
    end
end