function gierer_meinhardt
% MATLAB function to simulate the Gierer-Meinhardt reaction-diffusion model
% 

%% Model parameters
% diffusion coefficient
Du = 1.0;
Dv = 13;
Dv = 20;
% kinetic parameters
a = 0.2;
b = 3*a;

%% Simulation parameters
% domain length
length = 30;
% length = 5.3; % no patt
% number of grid points in the spatial discretisation
npts = 101;
% mesh spacing in the spatial discretisation
dx = length/(npts-1);
dx2 = dx^2;
% number of intermediate time plots
nplots = 30;
% final time
fin = 300;
% relative size of random spatial perturbation
mypert = 0.01;
% mypert = 0.5;
% mypert = 0.0;
%
% relative size of random homogeneous perturbations
upert =  0.;
vpert = 0.;

%% main code - should not need editing
% time between plots
tinterval = fin/nplots;
% steady states
us = a+b;
vs = b/us^2;

% trace and determinant of the Jacobian
tr = fu(us,vs)+gv(us,vs);
det = fu(us,vs)*gv(us,vs) - fv(us,vs)*gu(us,vs);
% for Turing instability we need this > 0
t1 = gv(us,vs) + Dv*fu(us,vs);
% and this < 0
t2 = 4*Dv*det - (gv(us,vs)+Dv*fu(us,vs))^2;

% Turing instability if Dv/Du > dcond
dcond = 32*a^2*(5+2*sqrt(6))

tcondstr = sprintf('Turing conditions.\n Need t1>0: t1=%g\n t2<0: t2=%g\n Equivalently need Dv-Du*dcond=%g>0',t1,t2,Dv-Du*dcond);
disp(tcondstr)

% a string to use in the figure title
figstr = 'Reaction diffusion simulation ';

% set the initial conditions for U and V
U0 = us*(ones(npts,1) + mypert*(2*rand(npts,1)-1) + upert*(2*rand-1));
V0 = vs*(ones(npts,1) + vpert*(2*rand-1));
ics = [U0; V0];

% set up the spatial mesh
xx = [0:dx:(npts-1)*dx];

% initialise time and the plot counter
i = 0;
t = 0;

% plot initial conditions
figure(1)
clf
subplot(2,1,1)
plot(xx,U0,'linewidth',2)
ylim([0 1.5])
set(gca,'fontsize',14)
title([figstr ': i = ',num2str(i),', time = ',num2str(t(end))])
hold on

subplot(2,1,2)
plot(xx,V0,'linewidth',2)
xlabel('x, distance','fontsize',16)
ylabel('V','fontsize',16)
ylim([0 1.1])
set(gca,'fontsize',14)
hold on

% wait a bit before continuing
pause(0.5)
%
U = ics(1:npts);
V = ics(npts+1:2*npts);

for i=1:nplots
    % plot the previous solution in blue
    subplot(2,1,1)
	plot(xx,U,'b','linewidth',2)
    ylabel('U','fontsize',16)
	title([figstr ': i = ',num2str(i),', time = ',num2str(t(end))])
    subplot(2,1,2)
	plot(xx,V,'b','linewidth',2)
    xlabel('x, distance','fontsize',16)
    ylabel('V','fontsize',16)

    tspan = [(i-1)*fin/nplots (2*i-1)*fin/(2*nplots) i*fin/nplots];
	options = [];
	[t,y] = ode23s(@meinhardt_ode,tspan,ics,options);
	ics = y(end,:)';
    U = ics(1:npts);
    V = ics(npts+1:2*npts);

    % plot the new solution in red - so the red curve is always the latest one
    subplot(2,1,1)
	plot(xx,U,'r','linewidth',2)
    ylabel('U','fontsize',16)
	title([figstr ': i = ',num2str(i),', time = ',num2str(t(end))])
    subplot(2,1,2)
	plot(xx,V,'r','linewidth',2)
    xlabel('x, distance','fontsize',16)
    ylabel('V','fontsize',16)
    
    % pause a bit between frames
    pause(0.1)
end

    % -----------------------------------------------------------------------
    % Nested function -- problem parameters provided by the outer function.
    %
    function out1 = meinhardt_ode(t,y,flag)

    % split the input vector y into its U and V components
    U = y(1:npts);
    V = y(npts+1:2*npts);

    % diffusion with 'no-flux' bcs
    diffU = (U([2,1:npts-1]) - 2*U + U([2:npts,npts-1]))*(Du/dx2);
    diffV = (V([2,1:npts-1]) - 2*V + V([2:npts,npts-1]))*(Dv/dx2);

    % specify dydt = [dUdt; dVdt] 
    out1 = [a - U + (U.^2).*V + diffU;
            b - (U.^2).*V + diffV];

    end

    % -----------------------------------------------------------------------
    % Entries in the Jacobian matrix. 
    %
    function f1 = fu(u,v)
        f1 = 2*u*v-1;
    end
    function f2 = fv(u,v)
        f2 = u^2;
    end
    function g1 = gu(u,v)
        g1 = -2*u*v;
    end
    function g2 = gv(u,v)
        g2 = -u^2;
    end

end