% The Mathematics of Patterns 2016
%  Schnakenberg kinetics in 1D
% Bivash Kaity
% This MATLAB code will show you how to numerically simulate the 1D 
% Schnakenberg kinetics system, and make the (simplified) version of the 
% animations in the Turing Instabilities II section of the website. 
clear all
close all

%%%%%%%%%%%%%%%%%%%%%
%%%    Set-up     %%%
%%%%%%%%%%%%%%%%%%%%%

% Parameter values
c1=0.1; c2=0.9; cm1=1.0; c3=1.0; 
Du = 1;   Dv = 40;  % Diffusion constants

% Grid and initial data:
% w = 10;  % no pattern 
w = 100;   % pattern

Nx = 1000; % How many points we want to discretize our domain with
x = linspace(0,w, Nx); 
dx = x(2) - x(1); 

dt = 0.1; % size of our time step
t = 0:dt:100;   
Nt = length(t); % Number of time points

% Set up for the surface 
[X, T] = meshgrid(x, t); 
U = 0*X;
V = 0*X;

% Easier to deal with column vectors
x = x(:);
t = t(:);

%Initial conditions: small perturbation away from steady state
%Steady state (u0,v0)=[(c1+c2)/cm1, c2*cm1^2/c3*(c1+c2)^2]=[1,0.9]
u = (c1+c2)*ones(length(x),1) + 0.01*rand(Nx, 1); 
v = (c2/(c1+c2)^2)*ones(length(x),1);

% Save initial conditions
U(1,:) = u;
V(1,:) = v;


%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Making the matrix %%%
%%%%%%%%%%%%%%%%%%%%%%%%%

% To begin, let us recall how to set up the matrices used in the explicit
% and implicit finite difference methods.


%%% Forward (explicit) method %%%
% We want a tridiagonal matrix (see notes for details) 
% % a = (1-2*Du*dt/dx^2);  % values along the diagonal
% % b = Du*dt/dx^2;        % values in the off-diagonal
% % main = a*sparse(ones(Nx,1)); 
% % off  = b*sparse(ones(Nx-1,1));
% % Bu = diag(main) + diag(off,1) + diag(off,-1); %Still a sparse matrix
% % % Satisfying boundary conditions (no flux boundary)
% % Bu(1,1) = (a+b);
% % Bu(end,end) = (a+b);
 
% To have a more numerically stable code, we use the implicit method. 

%%% Backward (implicit) method %%%
% For u
    a = (1+2*Du*dt/dx^2);  % values along the diagonal
    b = Du*dt/dx^2;        % values in the off-diagonal
    main = a*sparse(ones(Nx,1)); 
    off  = -b*sparse(ones(Nx-1,1));
    Bu = diag(main) + diag(off,1) + diag(off,-1); %Still a sparse matrix
%     Satisfying boundary conditions
%  
%     Bu(1,1) = a-b;
%     Bu(end,end) = a-b;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    Bu(1,2) = -(2*b);
    Bu(end,end-1) = -2*b;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     Bu(1,1) = a-b;
%     Bu(end,end-1) = -2*b;
% Same thing for v
    a = (1+2*Dv*dt/dx^2); b = Dv*dt/dx^2;
    main = a*sparse(ones(Nx,1));
    off  = -b*sparse(ones(Nx-1,1));
    Bv = diag(main) + diag(off,1) + diag(off,-1);
%     For boundary condition at x=0: u0=u1 & at x=L: u_(Nx+1)=u_(Nx-1)
%  i.e. at left boundary have taken difference between 0 & 1
% where as at right boundary taken central differnece (Nx+1) & (Nx-1)

%     Bv(1,1) = a-b;
%     Bv(end,end) = a-b;
    
    Bv(1,2) = -(2*b);
    Bv(end,end-1) = -2*b;
%     Bv(1,1) = a-b;
%     Bv(end,end-1) = -2*b;

%%%%%%%%%%%%%%%%%%%%%%%%%
%%%     Plotting      %%%
%%%%%%%%%%%%%%%%%%%%%%%%%

figure(1); %create new figure
plot(x,u,'g.-', 'linewidth',1);
hold on;
plot(x,v,'r.-', 'linewidth',1);
hold off;
%axis([-1 80 -.01 15.01])  % Fix axis limits

for j = 2:Nt
    % f and g are the reaction terms in the G-M system
    f = c1-cm1*u+(c3.*v.*u.^2);
    g = c2-(c3.*v.*u.^2);

    % At each step we need to solve the system
    u = Bu\(u + dt*f);   % backward Euler
    v = Bv\(v + dt*g);
    
    % Plot
    plot(x,u,'b.-', 'linewidth',1);
    hold on;
    plot(x,v,'r.-', 'linewidth',1);
    hold off;
    %axis([-1 80 -.01 15.01])
    title(['t = ', num2str(j*dt)],'fontsize',24)
    xlabel('x','FontSize', 30); ylabel('u/v','FontSize', 30)
    drawnow; 
    
    % Save for surface
    U(j,:) = u;
    V(j,:) = v;
end
%Plotting the surface %%%%
% figure(2);
% s = surf(x, t, U);
% set(s, 'EdgeColor', 'none', 'FaceColor', 'interp'); % Sets up the colors
% xlabel('x') 
% ylabel('t') 
% zlabel('u')

%contour plot %%%
figure(3);
p = pcolor(x, t, U);
xlabel('x','FontSize', 30); ylabel('t','FontSize', 30)
set(p, 'EdgeColor', 'none', 'FaceColor', 'interp');
