% mostly by cbm, bwo, 20060210
% u is used throughout instead of psi

global N L ep U LAM

k = 2;
l = 2;

ep = 1/16;  % epsilon

% there is a div by 4 here because our domain is [-1,1]x[-1,1]
% (to match Trefethen's code)
lam0 = -(k^2+l^2)*pi^2/4;

% chebychev stuff: from Trefethen's book via Todd
N = 20; [D,x] = cheb(N); y = x;
[xg,yg] = meshgrid(x(2:N),y(2:N));
xx = xg(:); yy = yg(:);
D2 = D^2; D2 = D2(2:N,2:N); I = eye(N-1);
L = kron(I,D2) + kron(D2,I);

% this is psi0
uu = 1*sin(k*(pi/2*xg + pi/2)) .* sin(l*(pi/2*yg + pi/2));
u = uu(:);

figure(3); clf;
surf(xg, yg, uu);
xlabel('x');
ylabel('y');
pause(0);

opt = optimset('Display', 'iter');
opt = optimset(opt, 'TolFun', 1e-12, 'TolX', 1e-12);
opt_lam = optimset(opt, 'MaxFunEvals', 1000, 'MaxIter', 1000);

uold = u;
lamold = lam0;

for j = 1:2  % number of iterations
  U = uold;   % pass the fixed u as a global var
  %[lamnew, fval, exf] = fsolve(@iter_f_lam, lamold, opt)
  [lamnew, fval, exf] = fminsearch(@iter_f_lam, lamold, opt_lam)
  difflam = abs(lamold-lamnew)

  % we use fsolve next because fminsearch doesn't to vectors
  LAM = lamnew;  % pass the fixed lambda as a global var
  [unew, fval, exf] = fsolve(@iter_f_psi, uold, opt);
  fvalnorm = norm(fval)
  udiffnorm = norm(unew-uold,2)/length(unew)
  udiffnorminf = max(abs(unew-uold))

  usol = reshape(unew, N-1, N-1);
  figure(4); clf;
  surf(xg, yg, usol)

  uold = unew;
  lamold = lamnew;
  disp('**************************************************');
  disp('press enter to finish this iteration');
  pause;
end

% check:
res = L*unew - lamnew*unew - ep*unew.^3;
norm(res,2)
norm(res,inf)
max(abs(res))