function meinhardtg6( nScenario)

%meinhardt: Simulate Meinhardt's (2001) n-th scenario.
%----------------------------------------------------------
% Here I have implemented Meinhardt's equations exactly as
% he shows them in his Basic program, in order to get you
% started on your project. However there are several things
% you must do, which I haven't done here:
% 1. Meinhardt's simulation is 1-D - you need 2-D.
% 2. This uses Euler integration of the DE's - you need RK.
% 3. Meinhan extremelrdt uses ay long time-step dt=1. You
%	 may well need a much smaller one, which will make
%	 your simulation slower.
% 4. Meinhardt uses Fixed boundary conditions - I have used
%	 Periodic boundary conditions. You may use either
%	 Periodic, Fixed or Zero-flux boundary conditions to
%	 achieve your target pattern. (See Miura & Maini (2004)
%	 for an explanation of these boundary conditions).
% 5. You will also need to match your pattern against a
%	 target pattern, and you will perform this comparison
%	 by using a Fourier classification of both your
%	 generated pattern and your target pattern.
% 6. Beware of ill-conditioning! You should be able to
%	 recognise it when it occurs, and hopefully using the
%	 RK method should help, but if not, you will need to
%	 reduce the length of your time-step to restore
%	 well-conditioned behaviour.
%----------------------------------------------------------
if nargin < 1
	nScenario = 1;	% Periodic pattern
end;

% Set up scenario constants:
switch nScenario
	case 1
		% Periodic pattern:
		Tmax = 70000;	Xmax = 80;			% Max t and x value
		dt = .1;		dx = 1;				% Step-sizes
		au = 0.010;	av = 0.061;              % Decay constants
		bu = 0.1;	bv = 0;                 % Base expression rates
		Du = 0.01;	Dv = 0.3;              % Diffusion constants
	otherwise
		disp( 'Sorry: unknown scenario!');
		return;
end;

% Define the domain:
Nt = floor(Tmax/dt) + 1;
Nx = floor(Xmax/dx) + 1;
T = dt*(0:(Nt-1));
X = dx*(0:(Nx-1));


% Set up slightly noisy initial conditions:
noise = 0.005;
lu = au*(1 + noise*(2*rand(Nx,Nx)-1));
u = ones(Nx,Nx);
v = u;

% Reaction rules:
f = @(u,v) lu.*(u.^2 + bu)./(v) - au*u;
g = @(u,v) lu.*u.^2 - av*v + bv;


% Display initial state:
figure(1);
pcoloru =  pcolor(X,X,u);
%pcoloru = surf(u);
shading interp;
colormap('jet');
figure(2);
pcolorv =   pcolor(X,X,v);
%pcolorv = surf(v);
shading interp;
colormap('jet');
drawnow;



% u(round(Nx/2),round(Nx/2)) = 100;
% %v((round(Nx/2)+15),round(Nx/2)) = 10000;
%  v(1:(round(Nx/2)),round(Nx/2)) = 1000;




% Simulation loop:
for t = 1:Nt-1

    u(round(Nx/2),round(Nx/2)) = 1;
    v(1:(round(Nx/2)),round(Nx/2)) = 5;
	u = u + dt * (f(u,v) + Du*coolshift2(u,dx));
    %v(round(5*end/1):round(6*end/11),end) = 20000;
    v = v + dt * (g(u,v) + Dv*coolshift2(v,dx));
    %v(1:(round(Nx/2)),round(Nx/2)) = 50;
    %v((end),end) = 20000;
   % v(round(4*end/9):round(5*end/9),end) = 2000;

    if mod(t,500) == 0
      t
    %  plot(X,u,X,v);
       set(pcoloru,'CData',u); % update with new data u
       set(pcolorv,'CData',v); % update with new data v
      drawnow;	% This is necessary if we want animation
    end


%       plot(X,u,'b', X,v,'r');
%       surf(u);
% 		shading flat;
%       if (mod(t,250) == 0)
%       drawnow;   % This is necessary if we want animation
%    end