%
%  loop control for lia's problem
%

n = 2*32;  h = 1/n;
A = 1;  B = exp(i*2*pi/3);   C = exp(i*4*pi/3); 

if (ei2 == 1)
  nonlin = zeros((n-2)^2,1);
end

disp(' ')
disp(['  ei2 = ' num2str(ei2)])
disp(' ')

%  first step (strip init line from my_delsqdemo)
G = numgrid('N',n);  N = sum(G(:)>0);  lapl = -delsq(G)/(h^2);

lia;  u_old = u;  pause(1)

%  iteration
change = 1;
while (change > 1e-8)
  nonlin = ( (abs((u-A).*(u-B)).^2).*(u-C) ...
           + (abs((u-B).*(u-C)).^2).*(u-A) ...
           + (abs((u-C).*(u-A)).^2).*(u-B) ) * ei2;

  lia;  change = abs(max(u-u_old))
  u_old = u;  pause(1)
end

%  plot
%    120 contour ghosts

figure(2);clf;  gray = 0.8*[1 1 1];
[nmin node]= min(diag(flipud(abs(U(2:end-1,2:end-1)))));
node = node+1;
trip1 = [node*(1+i)  (node*(1+i)+n/4*exp(i*(pi/4-2*pi/3)))];
trip2 = [node*(1+i)  (node*(1+i)+n/4*exp(i*(pi/4+2*pi/3)))];
plot(real(trip1),imag(trip1),'color',gray','linewidth',3);
hold on
plot(real(trip2),imag(trip2),'color',gray','linewidth',3);
plot(1:n,1:n,'color',gray','linewidth',3);

%  phase boundaries

contour(1:n,n:-1:1,angle(U*exp( i*pi/6))-pi/6,[-pi -pi/3 pi/3]);
axis([1 n 1 n]);  axis square;
caxis([-pi pi]);  colormap(hsv);  colorbar;

title('\bf phase boundaries')
xlabel('\bf x-axis matrix indices')
ylabel('\bf y-axis matrix indices')

figure(1)