% Test for solving the H-equation (Kelley book, p87)
% via L-Broyden's method
clear; close all
c = 1 - 1e-4; N = 4000;

x0 = ones(N,1);
tol = 5e-12;
maxiter = 100;
prt = 1;

run_my_codes = input('run_my_codes (default is false) = ');
if isempty(run_my_codes), run_my_codes = 0; end

if run_my_codes
    strB = 'myBroyden(@myHeqn,x0,tol,maxiter,prt,c,N);';
    strL = 'myLBroyden(@myHeqn,x0,tol,maxiter,L,prt,c,N);';
else % run TZ codes
    strB = 'yzBroyden(@yzHeqn,x0,tol,maxiter,prt,c,N);';
    strL = 'yzLBroyden(@yzHeqn,x0,tol,maxiter,L,prt,c,N);';
end

fprintf('Broyden: c = %g, N = %i\n',c,N)
t0 = tic;
[x,nrmF,iter] = eval(strB);
fprintf('Broyden: mean(x) = %.8f\n',mean(x))
fprintf('nrmF %g, iter %i, time %g\n\n',nrmF(end),iter,toc(t0))
h = semilogy(1:iter,nrmF,'-o'); set(h,'linewidth',2);

for L = 5:-1:3
    fprintf('L-Broyden: c = %g, N = %i  L = %i\n',c,N,L)
    t0 = tic; hold on;
    [x,nrmF,iter] = eval(strL);
    fprintf('L-Broyden using L = %i: mean(x) = %.8f\n',L,mean(x))
    fprintf('nrmF %g, iter %i, time %g\n\n',nrmF(end),iter,toc(t0))
    h = semilogy(1:iter,nrmF,'-o');
    set(h,'linewidth',2); grid on;
    title('Convergence History'); xlabel('Iter'); ylabel('F-norm');
end
legend('Broyden','L-Broy: L=5','L-Broy: L=4','L-Broy: L=3'); 
hold off; shg
