clear all;

men0 = 1;
hold_it = 0;
i_fig = 0;
cont_form = 1;
rlam = 1; slam = 1; principal_root = 1;  spurious_1  = 1; spurious_2 = 1; 
max_root = 1; title_meth = 1;
while men0 >= 0

  if men0 == 1
    meth = menu(...
	' Method choices ', ...		
	' Euler Explicit ', ...		% 1
	' Leapfrog ', ... 		% 2
	' AB2 ', ...			% 3
	' AB3 ', ...			% 4
	' Trapezoidal Implicit ', ...	% 5
	' Euler Implicit ', ...		% 6
	' 2nd Order Backward ', ...	% 7
	' AM3 ', ...			% 8
	' ABM3 ', ...			% 9
	' Gazdag ', ...		        %10
	' RK2 ', ...		        %11
	' RK4 ', ...		        %12
	' Milne 4th ', ...		%13
	' Midterm 1999 Question 3 ',...	%14
	' New figure', ...              %15
	' preferences',...              %16
	'exit');			%17
    if meth == 17 break; end
    if meth == 16
      mnn = menu('preferences','Contours','Contour = 1','toggle hold on');
      if mnn == 1 cont_form = 0;
      elseif mnn == 2 cont_form = 1;
      elseif mnn == 3 hold_it = 1-hold_it;
      end
    elseif meth == 15 
      i_fig = 0;  i_fig = figure;

    else
      % $$$   xsize = input('Enter max h lam ');
      xsize = 2.5;
      if meth == 2 | meth == 3 | meth == 4 | meth == 10 xsize = 1.5; end
      if meth == 6 xsize = 10; end
      if meth == 8 | meth == 7 xsize = 6.5; end
      if meth == 9 xsize = 4; end
      if meth == 12 xsize = 3; end
      if meth == 14 xsize = 1; end


      jmax = 200;
      clear rlam slam principal_root spurious_1 spurious_2 max_root title_meth;
      rlam = -xsize + 2.*xsize*[0:jmax-1]/jmax;
      slam = i*rlam;

      principal_root = zeros(jmax,jmax);
      spurious_1 = principal_root;
      spurious_2 = principal_root;
      max_root = principal_root;

      if meth == 1
	title_meth = 'Euler Explicit: One Root';
	for k = 1:jmax
	  for j = 1:jmax
	    hlam = rlam(k) + slam(j);
	    principal_root(k,j) = 1.0 + hlam;
	  end
	end
      elseif meth == 2
	title_meth = 'Leapfrog: Two Roots';
	for k = 1:jmax
	  for j = 1:jmax
	    hlam = rlam(k) + slam(j);
	    principal_root(k,j) = hlam + sqrt(1 + hlam^2);
	    spurious_1(k,j) = hlam - sqrt(1 + hlam^2);
	  end
	end

      elseif meth == 3
	title_meth = 'AB2: Two Roots';
	for k = 1:jmax
	  for j = 1:jmax
	    hlam = rlam(k) + slam(j);
	    principal_root(k,j) = 0.5*(1.0 + 3/2*hlam + sqrt(1 + hlam + 9/4*hlam^2));
	    spurious_1(k,j) = 0.5*(1.0 + 3/2*hlam - sqrt(1 + hlam + 9/4*hlam^2));
	  end
	end
      elseif meth == 4
	title_meth = 'AB3 : Three Roots';
	for k = 1:jmax
	  for j = 1:jmax
	    hlam = rlam(k) + slam(j);
	    rp4 = 1;
	    rp3 = -(1+23/12*hlam);
	    rp2 = 16/12*hlam;
	    rp1 = -5/12*hlam;
	    p = [rp4 rp3 rp2 rp1];
	    ss = roots(p);
	    size_s = size(ss);
	    principal_root(k,j) = max(abs(ss));
	  end
	end

      elseif meth == 5
	title_meth = 'Euler Implicit: One Root';
	for k = 1:jmax
	  for j = 1:jmax
	    hlam = rlam(k) + slam(j);
	    principal_root(k,j) = 1.0/(1.0-hlam);
	  end
	end

      elseif meth == 6
	title_meth = 'Trapezoidal Implicit: One Root';
	for k = 1:jmax
	  for j = 1:jmax
	    hlam = rlam(k) + slam(j);
	    principal_root(k,j) = (1.0+0.5*hlam)/(1.0-0.5*hlam);
	  end
	end
	
      elseif meth == 7
	title_meth = '2nd Order Backward: Two Roots';
	for k = 1:jmax
	  for j = 1:jmax
	    hlam = rlam(k) + slam(j);
	    principal_root(k,j) = (2/3 + 1/3*sqrt(1 + 2*hlam))/(1-2/3*hlam);
	    spurious_1(k,j) = (2/3 - 1/3*sqrt(1 + 2*hlam))/(1-2/3*hlam);
	  end
	end

      elseif meth == 8
	title_meth = 'AM3: Two Roots';
	for k = 1:jmax
	  for j = 1:jmax
	    hlam = rlam(k) + slam(j);
	    b_l = -(1.0 + 2/3*hlam);
	    a_l = 1.0 - 5/12*hlam;
	    c_l = 1/12*hlam;
	    principal_root(k,j) = (-b_l + sqrt(b_l^2 - 4.0*a_l*c_l))/(2.0*a_l);
	    spurious_1(k,j) = (-b_l - sqrt(b_l^2 - 4.0*a_l*c_l))/(2.0*a_l);
	  end
	end

      elseif meth == 9
	title_meth = 'ABM3: Two Roots';
	for k = 1:jmax
	  for j = 1:jmax
	    hlam = rlam(k) + slam(j);
	    b_l = -(1+13/12*hlam+15/24*hlam^2);
	    c_l = 1/12*hlam*(1+5/2*hlam);
	    a_l = 1;
	    principal_root(k,j) = (-b_l + sqrt(b_l^2 - 4.0*a_l*c_l))/(2.0*a_l);
	    spurious_1(k,j) = (-b_l - sqrt(b_l^2 - 4.0*a_l*c_l))/(2.0*a_l);
	  end
	end

      elseif meth == 10
	title_meth = 'Gazdag: Three Roots';
	for k = 1:jmax
	  for j = 1:jmax
	    hlam = rlam(k) + slam(j);
	    rp4 = 1;
	    rp3 = -(1+2*hlam);
	    rp2 = 3/2*hlam;
	    rp1 = -1/2*hlam;
	    p = [rp4 rp3 rp2 rp1];
	    ss = roots(p);
	    size_s = size(ss);
	    principal_root(k,j) = max(abs(ss));
	  end
	end

      elseif meth == 11
	title_meth = 'RK2: One Root';
	for k = 1:jmax
	  for j = 1:jmax
	    hlam = rlam(k) + slam(j);
	    principal_root(k,j) = 1 + hlam + 1/2*hlam^2;
	  end
	end

      elseif meth == 12
	title_meth = 'RK4: One Root';
	for k = 1:jmax
	  for j = 1:jmax
	    hlam = rlam(k) + slam(j);
	    principal_root(k,j) = 1 + hlam + 1/2*hlam^2 + 1/6*hlam^3 + 1/24*hlam^4;
	  end
	end
	
      elseif meth == 13
	title_meth = 'Milne 4th: Two Roots';
	for k = 1:jmax
	  for j = 1:jmax
	    hlam = rlam(k) + slam(j);
	    b_l = -4/3*hlam;
	    a_l = 1.0 -1/3*hlam;
	    c_l = -(1 + 1/3*hlam);
	    principal_root(k,j) = (-b_l + sqrt(b_l^2 - 4.0*a_l*c_l))/(2.0*a_l);
	    spurious_1(k,j) = (-b_l - sqrt(b_l^2 - 4.0*a_l*c_l))/(2.0*a_l);
	  end
	end

      elseif meth == 14
	title_meth = 'Midterm 1999 Question 3: Two Roots';
	for k = 1:jmax
	  for j = 1:jmax
	    hlam = rlam(k) + slam(j);
	    principal_root(k,j) = 0.5 + hlam + 0.5*sqrt(1 + 4.0*hlam^2);
	    spurious_1(k,j) = 0.5 + hlam - 0.5*sqrt(1 + 4.0*hlam^2);
	  end
	end

      end

      if meth < 15
% thp 	if i_fig == 0 & hold_it == 0
% thp 	  i_fig = figure;
% thp 	  set(i_fig,'Position',[300,100,800,800]);
% thp 	  clf;
% thp 	else
% thp 	  if i_fig == 0 i_fig = figure; else figure(i_fig); end
% thp 	end

        i_fig = figure;
	
	numcontours = 1;  minall_con = 1.0;
	
	x = zeros(jmax,jmax);
	y = x;
	x = rlam;
	y = imag(slam);
	z = max(abs(principal_root),abs(spurious_1));
	
	
	v = minall_con;
	V = v;
	c = contour(x',y',z',[v v],'k');
	hold on;
	if cont_form == 0
	  maxall_con = 2.0;
	  n_con = 6;
	  delcontours = (maxall_con-1.00)/(n_con-1);
	  v = 1.00 + [0:n_con-1]*delcontours;
	  v = round(100.*v)*0.01;
	  c = contourf(x',y',z',v);
	end
	% $$$     grid on;
	axis('square');
	title(title_meth);
	xlabel('Real(h \lambda)');
	ylabel('Imag(h \lambda)');
	
      end
    
    end  
  
  end

  men0 = menu('Stability in Lambda Plane','Another method','Add Lambda spectrum','Exit');
  if men0 == 3 break; end
  
  if men0 == 2
    hold on;
    Lambda_spectrum(i_fig);
    men0 = 1;
  end
% thp   if men0 == 1
% thp     clear rlam slam principal_root spurious_1 spurious_2 max_root title_meth;
% thp     men0 = 1;
% thp   end
  

end