function lambda_plane_eval_search

clear; clf;
L=10;           %U-integration start
xspan=[L -L];  %x integration range
numsteps=100;  %number of steps used in calculating U_lambda

%setup graphics boxes
eigenvalues=[-0.75,0,1.25];
subplot(211);
plot3(real(eigenvalues), imag(eigenvalues), zeros(3),'.k');
hold on; grid on; title('\lambda paths and centres in C-plane');
xlabel Re(\lambda); ylabel Im(\lambda); zlabel '\Delta \theta(s)/(2\pi)';
subplot(212); title('\theta(s)/(2\pi) for the various paths'); hold on; grid on; ylabel \theta(s)/(2\pi); xlabel s;
f=getframe;

lambda_0_RE_counter=0;

for lambda_0_RE=-0.80:0.20:1.6
    lambda_0_RE_counter = lambda_0_RE_counter+1;
    lambda_0_IM_counter=0;
    
    for lambda_0_IM = -0.4:0.20:0.4
        lambda_0_IM_counter = lambda_0_IM_counter+1;

        %lambda path definition along loop around lambda_0
        lambda_0 = lambda_0_RE + i*lambda_0_IM; %centre of lambda path circle 
        r=0.12;                                  %radius of lambda path circle

        subplot(211)
        for j=1:numsteps
            lambda_array(j)=lambda_0 + r*exp(2*pi*i*j/numsteps);
        end
        plot3(real(lambda_array),imag(lambda_array),zeros(length(lambda_array)),'r');
        plot3(real(lambda_0), imag(lambda_0),0,'.r');
        %f=getframe;

        %lambda path integration
        theta_init=0;           %theta initial value
        sspan=[0 2*pi];
        [s theta] = ode45(@theta_integration, sspan, theta_init);

        subplot(212)
        plot(s,theta/(2*pi))
        delta_theta = (theta(length(theta))-theta(1))/(2*pi);
        subplot(211)
        plot3(real(lambda_0),imag(lambda_0),delta_theta,'*k')
        Z(lambda_0_IM_counter,lambda_0_RE_counter)=delta_theta;
    end
end

figure
meshc(Z)
xlabel 'Re(\lambda) (-1->2)'; ylabel 'Im(\lambda) (-1->1)'; zlabel '\Delta \theta(s)/(2\pi)';


%Nested Functions
%-----------------------------------------------------------
    function theta_s = theta_integration(s,theta)
        
        next_lambda=lambda_0 + r*exp(i*(s+(2*pi/numsteps)));
        present_lambda = lambda_0 + r*exp(i*s);
        last_lambda=lambda_0 + r*exp(i*(s-(2*pi/numsteps)));

        next_mu=sqrt(next_lambda+1);
        present_mu=sqrt(present_lambda+1);
        last_mu=sqrt(last_lambda+1);

        %plot current lambda
%        subplot(211)
%        plot3(real(present_lambda), imag(present_lambda),0,'.')
%        f=getframe;
                
        %next_lambda... integration of u from x=L -> 0
        lambda=next_lambda; mu=next_mu; ksi_1=[ 1; -mu]; %stable evec of A_infinity exp(-mu*L)
        u_init=ksi_1; extended_u_init = [u_init; theta_init]; %starting vector (u(L,lambda), theta_init)
        [x extended_u] = ode45(@u_ODE, xspan, extended_u_init);
        u=[extended_u(:,1) extended_u(:,2)];
        theta=[extended_u(:,3)]; next_u=u(length(u),:).'; %results of integration
        clear u extended_u

        %present_lambda... integration of u from x=L -> 0
        lambda=present_lambda; mu=present_mu; ksi_1=[ 1; -mu]; %stable evec of A_infinity exp(-mu*L)
        u_init=ksi_1; extended_u_init = [u_init; theta_init]; %starting vector (u(L,lambda), theta_init)
        [x extended_u] = ode45(@u_ODE, xspan, extended_u_init);
        u=[extended_u(:,1) extended_u(:,2)];
        theta=[extended_u(:,3)]; present_u=u(length(u),:).'; %results of integration
        clear u extended_u
        
        %last_lambda... integration of u from x=L -> 0
        lambda=last_lambda; mu=last_mu; ksi_1=[ 1; -mu]; %stable evec of A_infinity exp(-mu*L)
        u_init=ksi_1; extended_u_init = [u_init; theta_init]; %starting vector (u(L,lambda), theta_init)
        [x extended_u] = ode45(@u_ODE, xspan, extended_u_init);
        u=[extended_u(:,1) extended_u(:,2)];
        theta=[extended_u(:,3)]; last_u=u(length(u),:).'; %results of integration
        clear u extended_u
        
        u_lambda = (next_u - last_u)/(next_lambda - last_lambda);
        theta_s = (imag((present_u')*(u_lambda)*(i*r*exp(i*s))))/(present_u'*present_u);

        %----------------------------------------------------------
        function extended_u_x = u_ODE(x,extended_u)
            u = [extended_u(1)      %ie the real u column vector
                 extended_u(2)];
            f = (3/2)*((sech(x/2))^2);
            A = [	0                    1
                    lambda+1-2*f    0	];
            dtheta_dx = (imag(u'*A*u))/(u'*u);
            u_x = (A - (-mu)*eye(2))*u;
            extended_u_x = [u_x
                            dtheta_dx];
        end
        %----------------------------------------------------------
        
    end

%----------------------------------------------------------
end