% This file solves the diffusion equation with exponential growth or decay

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

cmykflg=1;
%if cmykflg=1, this code will make cmyk figure files
%cmykflg = 1;



L = 40; % length of the domain
d = 1; % diffusion coefficient
Kvals = [0.2,0,-0.2]; % values of k for this simulation

N=200; % number of spatial grid points
M = 12; %number of time step
h = L/N; %spatial discretization
dt = h %time step
rat = ones(M,3);
 scal = d*dt/h^2;

X = h*(1:N)'; % x grid coordinates
 
kstep = 6; 
kstep*dt%intermediate steps to take
t=0; %initialize time
tn = zeros(M,1);
 
% set up the CN matrix

Atm=(1+scal)*ones(N,1);
Atm(1,1) = 1+scal/2;
Atm(N,1) = 1+scal/2;

Am=diag(Atm)-diag(scal/2*ones(N-1,1),1)-diag(scal/2*ones(N-1,1),-1);

  V = 0.2*exp(-(X-L/2).^2/(L/4));
  U = V*ones(1,3);
   
% uses Crank Nicolson to solve the diffusion equation

for j = 1:3
    rat(1,j) = U(N/2,j)/U(N/2,2); % the ratio of two solutions
    figure(j)
  plot(X,U(:,j),'--','linewidth',2)
  axis([0 L  0 1])
  hold on
end

% The algorithm
 

for n=2:M
for k = 1:kstep
    
    for j = 1:3
        K = Kvals(j); % do all 3 simulations simultaneously
        V = U(:,j);
        
%  F = K*V.*(1-V).*(V-0.1);bistable equation;
% F = K*V.*(1-V); %Fisher equation
 F =  K*V; % simple growth
V_t =  scal/2 * ([V(2:N );V(N )] - 2*V + [V(1 );V(1:N-1 )]) +dt*F ;
U(:,j)  = Am\(V+V_t);
t=t + dt;

    end
 

end
for j = 1:3
    rat(n,j) = U(N/2,j)/U(N/2,2);
    figure(j)
  plot(X,U(:,j),'linewidth',2)
  axis([0 L  0 1])
 % axis([0 L -.5 1])
 xlabel('x', 'fontsize',20)
 ylabel('U(x,t)','fontsize',20)
 
hold on
end
tn(n) = t;

end

 if(cmykflg==1)
     figure(3)
print('../../figs_c/chapt_6/diff_w_dk','-deps','-cmyk')
figure(1)
print('../../figs_c/chapt_6/diff_w_grth','-deps','-cmyk')

  end

   
 figure(4)
 semilogy(tn,rat,'linewidth',2)
 text(25,4,'\sigma \alpha=0.2','fontsize',20)
 text(25,0.25,'\sigma \alpha=-0.2','fontsize',20)
 xlabel('t', 'fontsize',20)
 ylabel('ln(U(x,t)/u_0(x,t)','fontsize',20)
TL = t;
if(cmykflg==1) 
print('../../figs_c/chapt_6/diff_rts','-deps','-cmyk')
end
     