% Turing instability in Schnackenberg system
% 
% u_t =     u_xx + alpha*(a - u + u^2*v)
% v_t = d * v_xx + alpha*(b- u^2*v)
%       f(u,v) = (a - u + u^2*v), g(u,v) = (b- u^2*v).
% for 0 < x < 1, with u_x=v_x=0 at x=0,1.
% with d > 0, alpha > 0.
%
% solve by splitting --
%  first
%     u_t = u_xx
%     v_t = d c_xx
%
% then
%     u_t =  alpha*(a - u + u^2*v)
%     v_t =  alpha*(b- u^2*v)
% 
% 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.
%
% We want system to be stable to spatially uniform
% perturbations so we want tr(J) < 0 and det(J) >0.
% Since tr(J)= ((b-a) - (b+a)^3)/(b+a) and
%      det(J)= (b+a)^2,
% these conditions hold if (b-a) < (b+a)^3;
% 
% We require b > a so that f_u > 0.
%
% Use discrete approximations with spatial points xj = (j-1)*dx for
% j=1,...,N+2  where (N+1)*dx=L, i.e. (N+1)=L/dx
%
% We want to find approximations to the solution of the PDE
% at the N+2 points xj for j=1,...,N+2;
%
% The linear systems are solved with a direct method (using \) if necessary.
%
% define parameters
%
% domain length
  L = 1.0;
%  
  d    = input(' enter diffusion ratio d ');
  a    = input(' enter parameter a ');
  b    = input(' enter parameter b ');
  alpha= input(' enter parameter alpha ');
  N    = input(' enter number of interior spatial grid points N ');

  tend = input(' enter final time tend ');
  plotfreq = input(' time interval between plots ');
  lambda = input(' enter ratio lambda ');  
  random = input(' type 1 for random perturbation ');
  tplot  = min(plotfreq,tend);
%
% DEFINE NUMERICAL PARAMETERS
%
  dx = L/(N+1);
  dt = lambda*(dx*dx)/max(1,d); 
  ntmax = ceil(tend/dt)+1;
  lambdau =   dt/(dx*dx);
  lambdav = d*dt/(dx*dx);  
  x    = dx*[0:1:N+1]';
%
% SET UP SOLUTION VARIABLE ARRAYS
%
  u_old = zeros(N+2,1);
  u_new = zeros(N+2,1);
  v_old = zeros(N+2,1);
  v_new = zeros(N+2,1);
%  
  LEFT = 1;
  RIGHT = N+2;  
  dim=0;  
%
% DEFINE MATRICES NEEDED FOR CN METHOD
%
  [MLu,MRu] = matrixdefineTuring(N,lambdau,x,dx,dim);  
  [MLv,MRv] = matrixdefineTuring(N,lambdav,x,dx,dim);  
%
% DEFINE INITIAL CONDITIONS
%
  perturbsize = 0.1
  if random == 1,
     rand('state',sum(100*clock));
     upert(LEFT:RIGHT) = perturbsize*(rand(RIGHT-LEFT+1,1)-0.5);
     vpert(LEFT:RIGHT) = perturbsize*(rand(RIGHT-LEFT+1,1)-0.5);
  else
     xtemp = x(LEFT:RIGHT);
     upert(LEFT:RIGHT) = perturbsize*cos(pi*xtemp/L)+ ...
	 perturbsize*cos(2*pi*xtemp/L)+...
	 perturbsize*cos(3*pi*xtemp/L) + perturbsize*cos(4*pi*xtemp/L)+...
	 perturbsize*cos(5*pi*xtemp/L) + perturbsize*cos(6*pi*xtemp/L)+...
	 perturbsize*cos(7*pi*xtemp/L) + perturbsize*cos(8*pi*xtemp/L)
     vpert(LEFT:RIGHT) = 0.0;
  end;
%  
  ubar = (a+b);
  vbar =  b/(a+b)^2;  
  u_init(LEFT:RIGHT) = ubar + upert(LEFT:RIGHT);
  v_init(LEFT:RIGHT) = vbar + vpert(LEFT:RIGHT);
  u_old (LEFT:RIGHT) = u_init(LEFT:RIGHT);
  v_old (LEFT:RIGHT) = v_init(LEFT:RIGHT);
%  
% START FIGURE
%
  figure
  jj = 1
  subplot(5,2,jj)
  bar = max(ubar,vbar)
  plot(x,u_old,'b-'),axis([0 L 0.00 2.5*bar])
  hold on;
  plot(x,v_old,'r-')
  drawnow
  tic;
  tstart=cputime;
%
% MARCH FORWARD IN TIME
%
  t = dt;
  while(t < tend);
    umax = max(u_old(:)); vmax = max(v_old(:));    
    [t,umax,vmax]
%
% DO DIFFUSION FIRST
%
    bcvector = zeros(N+2,1);
    [u_new(LEFT:RIGHT)]=matsolveD(MLu,MRu,u_old(LEFT:RIGHT),bcvector);
    bcvector = zeros(N+2,1);
    [v_new(LEFT:RIGHT)]=matsolveD(MLv,MRv,v_old(LEFT:RIGHT),bcvector);
%
% DO REACTION SECOND -- 
%    
% Use 4th order Runge-Kutta ODE method
%
% For dy/dt = F(y,t), in each time step, do 
%      k1 = F(yn      ,tn   )
%      k2 = F(yn+.5*dt*k1,tn+.5*dt)
%      k3 = F(yn+.5*dt*k2,tn+.5*dt)
%      k4 = F(yn+   dt*k3,tn+   dt)
%      ynp1 = yn + dt/6 *(k1 + 2*k2 + 2*k3 + k4)
%
% Stage 1    
       utemp = u_new(LEFT:RIGHT);
       vtemp = v_new(LEFT:RIGHT);
       [k11(LEFT:RIGHT),k12(LEFT:RIGHT)] = reactionsTuring(utemp,vtemp,a,b,alpha);
% Stage 2
       utemp = u_new(LEFT:RIGHT) + .5*dt*k11';
       vtemp = v_new(LEFT:RIGHT) + .5*dt*k12';       
       [k21(LEFT:RIGHT),k22(LEFT:RIGHT)] = reactionsTuring(utemp,vtemp,a,b,alpha);
% Stage 3      
       utemp = u_new(LEFT:RIGHT) + .5*dt*k21';
       vtemp = v_new(LEFT:RIGHT) + .5*dt*k22';              
       [k31(LEFT:RIGHT),k32(LEFT:RIGHT)] = reactionsTuring(utemp,vtemp,a,b,alpha); 
% Stage 4      
       utemp = u_new(LEFT:RIGHT) + dt*k31';
       vtemp = v_new(LEFT:RIGHT) + dt*k32';                     
       [k41(LEFT:RIGHT),k42(LEFT:RIGHT)] = reactionsTuring(utemp,vtemp,a,b,alpha);
% Combine results for stages 1 to 4
       u_new(LEFT:RIGHT) = u_new(LEFT:RIGHT) + dt/6 * (k11 + 2*k21 + 2*k31 + k41)';
       v_new(LEFT:RIGHT) = v_new(LEFT:RIGHT) + dt/6 * (k12 + 2*k22 + 2*k32 + k42)';
       if(t > tplot)
	 jj = jj + 1
         tplot = tplot+plotfreq
         subplot(5,2,jj)	 
         plot(x,u_new,'b-'),axis([0 L 0.00 2.5*bar])
	 hold on
	 plot(x,v_new,'r-')
	 drawnow
      end;
      u_old = u_new;
      v_old = v_new;
      t = t + dt;
  end;
  toc
  tfinish=cputime;
  processor_time=tfinish-tstart
  name=strcat(' Schnackenberg a= ',num2str(a),' b= ',num2str(b),' d=',num2str(d),' \alpha=',num2str(alpha))
%  title(name)
%  xlabel('x')
%  ylabel('u,v')