function barrier3
global t
t = 1;
mu = 2;
x = [-3;3];
while (1/t > 1e-6)
    [x, fval, exitflag] = fsolve(@gf, x, optimset('Display','off'));
    fprintf(1,'t = %d, solution x = %8.6f, y = %8.6f\n',t, x(1), x(2));
    if (exitflag~=1)&(exitflag~=2)
        disp('Nonlinear solver does not converge! Choose a smaller mu.');
        return;
    end
    t = mu*t;
end


function [vv] = gf(vx)
global t
x = vx(1);
y = vx(2);
er = exp( - x^2 - y^2);
g = 2 - ( 0.5*x*y + (x+2)^2 + 0.5*(y-2)^2 ) ;
vv = [t*( er - 2*x^2*er + 0.1*x ) - 1/g * (-0.5*y-2*x-4); ...
    t*( -2*x*y*er + 0.1*y ) - 1/g * (-0.5*x-y+2) ];
t = 1, solution x = -2.091411, y = 2.138634
t = 2, solution x = -1.771123, y = 1.660566
t = 4, solution x = -1.282285, y = 0.922205
t = 8, solution x = -1.109355, y = 0.666282
t = 16, solution x = -1.041626, y = 0.568136
t = 32, solution x = -1.007943, y = 0.519471
t = 64, solution x = -0.990637, y = 0.494459
t = 128, solution x = -0.981780, y = 0.481649
t = 256, solution x = -0.977286, y = 0.475147
t = 512, solution x = -0.975021, y = 0.471868
t = 1024, solution x = -0.973883, y = 0.470221
t = 2048, solution x = -0.973313, y = 0.469396
t = 4096, solution x = -0.973028, y = 0.468983
t = 8192, solution x = -0.972885, y = 0.468776
t = 16384, solution x = -0.972814, y = 0.468673
t = 32768, solution x = -0.972778, y = 0.468621
t = 65536, solution x = -0.972760, y = 0.468595
t = 131072, solution x = -0.972751, y = 0.468582
t = 262144, solution x = -0.972747, y = 0.468576
t = 524288, solution x = -0.972744, y = 0.468573