%
%  solving the lia equation   del^2 u = C |(z-A)(z-B)|^2 (z-C) + ...
%    - with u given on bndry as in the assignment!
%    - a hacked version of "delsqdemo.m"
%

clf reset

% solve the Laplace boundary value problem:
%    delsq(u) = nonlin in the interior
%    u = A,B,C on the boundary

rhs = nonlin;
rhs(round(G(  2,2:n-1))) = rhs(round(G(  2,2:n-1))) -B/(h^2);
rhs(round(G(n-1,2:n-1))) = rhs(round(G(n-1,2:n-1))) -C/(h^2);
rhs(round(G(2:n-1,  2))) = rhs(round(G(2:n-1,  2))) -C/(h^2);
rhs(round(G(2:n-1,n-1))) = rhs(round(G(2:n-1,n-1))) -A/(h^2);

%  do the solve

flops(0);  tic; 
u = lapl\rhs;
disp(' ')
disp(['  # of flops = ' num2str(flops) '   ;   cpu time = ' num2str(toc)])
disp(' ')

%  plot the solution (with BCs & corners = 0)

U = complex(G);  
U(G>0) = full(u(G(G>0)));
U(1,2:end-1) = B;  U(2:end-1,1) = C;
U(n,2:end-1) = C;  U(2:end-1,n) = A;

%image(real(U),'cdatamapping','scaled');  colormap(jet);  
%image(imag(U),'cdatamapping','scaled');  colormap(jet);  
%image( abs(U),'cdatamapping','scaled');  colormap(jet);  
image(angle(U),'cdatamapping','scaled');  colormap(hsv);

colorbar;  axis square;  axis image
title('\bf complex triple point problem')
xlabel('\bf x-axis matrix indices')
ylabel('\bf y-axis matrix indices')