%
%     u_t =  alpha*(a - u + u^2*v)
%     v_t =  alpha*(b- u^2*v)
% 
  a    = input(' enter parameter a ');
  b    = input(' enter parameter b ');  
%
% uniform steady state is (ubar, vbar) where
%
  ubar = (a+b);
  vbar =  b/(a+b)^2;
%
% Require b>0 and (a+b) > 0
%
% f_u = (b-a)/(b+a); f_v = (a+b)^2;
% g_u = -2b/(a+b)  ; g_v = -(a+b)^2.
%
fu = (b-a)/(b+a)
fv = (a+b)^2
gu = -2*b/(a+b)
gv = -(a+b)^2
fu = -1 + 2*ubar*vbar
fv = ubar^2
gu = -2*ubar*vbar
gv = -ubar^2
trJ  = fu+gv
detJ = fu*gv - fv*gu
%
% The equation gamma_min=0 determines
% equation for critical diffusion coefficient
% fu^2 d^2 + (2*fu*gv - 4*detJ)*d + gv^2 = 0.
%
A = fu^2;
B = (2*fu*gv - 4*detJ);
C = gv^2;
%
% gamma_min < 0 for d > dp
%
dp = ( -B + sqrt(B^2 - 4*A*C))/(2*A)
dm = ( -B - sqrt(B^2 - 4*A*C))/(2*A)
%
% Check stability of mode cos(q*x); s = q^2;
flag = 1;
pi2 = pi*pi;
gammin = -5;
sigmax = 5;

nmax = input(' enter max mode number n to consider ')
figure
figure
while flag == 1
  d = input(' enter value of d ')
  alpha = input(' enter value of alpha ')
  n = [0:0.01:nmax];
  s = n.^2*pi2;  
  figure(1)
  plot([0 n(end)],[0 0],'r-.')
  hold on  
  figure(2)
  plot([0 n(end)],[0 0],'r-.')
  hold on
  [fu,fv,gu,gv];  
  [sig,bet,y] = schnackgamsig(s,fu,fv,gu,gv,alpha,d);
%  
% sig holds variable that we called sigma+ in class
% y   holds variable that we called gamma in class
%
  figure(1)
  plot(n,sig,'b-')
  sigmaxnew = max(sig);
  sigmax = max(sigmax,sigmaxnew);
  axis([0 n(end) -5 sigmax])
%%%  axis([0 n(end) -5 5])  
  hold on
  figure(2)
  plot(n,y,'r-')
  gamminnew = min(y);
  gammin = min(gammin,gamminnew);  
  if gammin < -50
     gammin = -50
  end
  axis([0 n(end) gammin 5])
%%%  axis([0 n(end) -250 50])  
  flag = input(' type 1 to continue ')
end;
name = strcat(' Schnackenberg a= ',num2str(a),' b= ',num2str(b),' d=',num2str(d))
%%,' \alpha = ',num2str(alpha))

name = strcat(' Schnackenberg a= ',num2str(a),' b= ',num2str(b),' \alpha = ',num2str(alpha))
figure(1)
title(name)
xlabel('n','Fontsize',16)
ylabel(' \sigma(q^2) -- q^2 = (n\pi)^2','Fontsize',16)
figure(2)
title(name)
xlabel('n','Fontsize',16)
ylabel(' \gamma(q^2) -- q^2 = (n\pi)^2','Fontsize',16)