Penalizacion externa

[x,f] = forPenExt([0,0],0.1,2,11,1*10^-6);
function f = penExt(x, r, q)
f = sin(x(2))*exp((1-cos(x(1)))^2) + cos(x(1))*exp((1-sin(x(2)))^2) + (x(1)-x(2))^2 ...
+r*((max(0,(x(1)+5).^2+(x(2)+5).^2-25)).^q);
end

function f = Func(x)
f = sin(x(2))*exp(1-cos(x(1)))^2 + cos(x(1))*exp(1-sin(x(2)))^2;
end

function [X,f] = forPenExt(X,r,q,c,e)
    index = 1;
    [x_int, fx_int] = fminunc(@(x)penExt(x,r,q),X);
    disp("****************************************");
    disp(index);
    disp(X);
    disp(x_int);
    disp(r);
    disp(fx_int);
    rk_int = r;
    x_old_int = X;
    disp(norm(abs(x_int - x_old_int)))
    disp("****************************************");
    while(e < norm(abs(x_int - x_old_int)))
        x_old_int = x_int;
        rk_int = rk_int*c;
        disp(index);
        disp(x_old_int);
        [x_int, fx_int] = fminunc(@(x)penExt(x,rk_int,q),x_old_int);
        disp(x_int);
        disp(rk_int)
        disp(fx_int);
        disp(norm(abs(x_int - x_old_int)))
        disp("****************************************");
        index = index+1;
    end
    X = x_int;
    disp(Func(X));
    f = Func(X);
end