% This file solves the system of diffusion-advection-reaction 
% equations for chemotaxis

clear
 
set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.2, ...
   'defaultlinelinewidth', 2.0, ...
   'defaultpatchlinewidth', 0.7);

 

L = 10; % length of the domain
d = 0.25; % diffusion coefficient of cAMP
mu = 0.07; % diffusion coefficient for dd
rho = 0.2; % rate of cAMP production
K = 0.3; % rate of cAMP degradation
chi = 0.2; % chemotactic parameter

N=201;    % number of spatial grid points
M = 12; % number of time steps
h = L/(N-1); %grid size

X = h*(1:N)' ;% grid points
 
 
dt = h;  % set time step
scald = mu*dt/h^2;
scalc = d*dt/h^2;
ast = 1;

 U = ast*ones(N,1)+ 0.1./cosh(10*(X-L/3.9)) ;% specify initial data for dd 
C = rho*U/K;% initial cAMP concentration
figure(1)
  plot(X,U,'--', 'linewidth',2 )
  
  hold on
kstep = 2000;
t=0;

% uses Crank Nicolson to solve the diffusion equation

% set up the CN matrices for C and U

Atmc=(1+scalc)*ones(N,1);
Atmc(1,1) = 1+scalc/2;
Atmc(N,1) = 1+scalc/2;

Amc=diag(Atmc)-diag(scalc/2*ones(N-1,1),1)-diag(scalc/2*ones(N-1,1),-1);


Atmd=(1+scald)*ones(N,1);
Atmd(1,1) = 1+scald/2;
Atmd(N,1) = 1+scald/2;

Amd=diag(Atmd)-diag(scald/2*ones(N-1,1),1)-diag(scald/2*ones(N-1,1),-1);

 
% The algorithm
 
for n=2:M
for k = 1:kstep
    % here is the code that does upwinding:
vel = chi*(C(2:N)-C(1:N-1))/h;  %calculate the velocity  % C is the chemoattractant
msk = (vel>0);
 Flx = [0;vel.*(msk.*U(1:N-1)+(1-msk).*U(2:N));0];% calculate the upwinded flux
Cdiff = (Flx(2:N+1)-Flx(1:N))/h;  % this is the gradient of the flux
U_t =  scald/2 * ([U(2:N );U(N)] - 2*U + [U(1);U(1:N-1)])-dt*Cdiff ;
C_t =  scalc/2 * ([C(2:N );C(N)] - 2*C + [C(1);C(1:N-1)])+dt*(rho*U-K*C) ;
  
U = Amd\(U+U_t);
C = Amc\(C+C_t);

t=t + dt;
end

  %axis([0 L  0 10])
 % axis([0 L -.5 1])
 figure(1)
  plot(X,U ,  'linewidth',2 )
   
xlabel('x','fontsize',20)
ylabel('u(x,t)','fontsize',20)
axis([0 L  0 18]) 

  
     
%legend('u','c','fontsize',20)
figure(2)
  plot(X,C ,  'linewidth',2 )
   
xlabel('x','fontsize',20)
ylabel('c(x,t)','fontsize',20)
axis([0 L  0 18]) 
 hold on
 
  
end
hold off

