% Brusselator (Ex.4 Kuznetsov 2nd ed., p. 32)
% Ver. 2 (selecting a and b with mouse)
  close all
  clear all

% Particular cases:
%   a = 2 ; b = 0.25 ; % ZONE I 
%   a = 1 ; b = 1.75 ; % ZONE II
%   a = 1 ; b = 2.5 ; % ZONE III (unstable)

%  Open window to select a and b with mouse
   figure (2)
   x1 = [0:.01:2] ; y1 = x1.^2 +1 ; y2 = (x1-1).^2 ;
   plot(x1,y1,'-k','linewidth',1.25) ; drawnow ; hold on
   plot(x1,y2,'-k','linewidth',1.25) ; drawnow ; 
   xlabel('a','fontsize',14)
   ylabel('b','fontsize',14)
   title('Please left-click on a point of the (a,b)-plane')
   axis equal
   axis ([0 2 0 2])
   [a,b] = ginput(1); 
   plot(a,b,'ok','markerfacecolor',[.5 .5 .5])
   text(a*1.02,b*1.02,['a = ', num2str(a), ',  b = ', num2str(b)])


% Time integration specs: dt = time-step ; nstep = # steps
  dt = 5e-2 ; nstep = 1500 ;

   figure(1)

% Grid of initial conditions
  for x0 = [0:1.5:3]
   for y0 =  [0:1.5:3]

   t = 0 ; x = [x0 ; y0] ; xvec = [x0] ; yvec = [y0] ; tvec = [t] ; 
  
  for tstep = 1:nstep
% Forward step via RK4
  p1 = brusfun(x,a,b) ;
  p2 = brusfun(x+0.5*dt*p1,a,b) ;
  p3 = brusfun(x+0.5*dt*p2,a,b) ;
  p4 = brusfun(x+dt*p3,a,b) ;
  
  x = x + (dt/6)*(p1+2*p2+2*p3+p4) ;
  t = dt*tstep ;
  xvec = [xvec x(1)] ; yvec = [yvec x(2)] ; tvec = [tvec t];
  end

  
  plot(xvec,yvec,'-k','linewidth',1.25) ; drawnow ; hold on
  plot(xvec(1),yvec(1),'ok','markerfacecolor',[1 1 1]) 

end
end

% Plotting steady sol. in phase plane
  plot(a,b/a,'ok','markerfacecolor',[0 1 0])
  text(a*1.01,(b/a)*1.01,'Equilibrium solution') 


 axis equal
 axis([0 4 0 4]) ; 
   
