%LogisticReactionDiffusion

n=100; % Number of spatial points
m=4000; % Number of temporal point % Final time point
D=.001; % Diffusion constant

%Construcstion of spatial and temporal domains and spacing.
x=linspace(0,1,n);
dx=x(2)-x(1);
dt=dx^2/(4*D); %This satisfies the CFL condition
t=linspace(0,(m-1)*dt,m);

%Construction of second derivative operator
DXX=-2*diag(ones(1,n),0)+diag(ones(1,n-1),1)+diag(ones(1,n-1),-1);
DXX(1,:)=0;
DXX(n,:)=0;
DXX=DXX./dx^2;

DXX=sparse(DXX);

% Construction of initial distribution
u0=sin(pi*x).*(sin(5*pi*x).^2);
v0=sin(pi/2).*(sin(2*pi*x).^2);
u0=u0';
v0=v0';


%% Method of lines
u=zeros(n,m); %Initialization of distribution
u(:,1)=u0;

v=zeros(n,m); %Initialization of distribution
v(:,1)=v0;


%Plotting Initial guess
figure('units', 'normalized', 'outerposition', [0,0,1,1]);
hold on;
plot(x,u0,x,v0,'LineWidth',2);
ylabel('$u(x,t)$','Interpreter','latex','Fontsize',18);
xlabel('$x$','Interpreter','latex','Fontsize',18);
title(strcat('$t=$',num2str(t(1))),'Interpreter','latex','Fontsize',18 );
axis([0,1,0,1]);
drawnow;

% Looping through Changes in time
for i=1:m-1
    u(:,i+1)=u(:,i)+(D*DXX*u(:,i)+u(:,i).*(1-u(:,i))-u(:,i).*v(:,i))*dt;
    v(:,i+1)=v(:,i)+(D*DXX*v(:,i)+v(:,i).*(1-v(:,i))-u(:,i).*v(:,i))*dt;
    
    u(1,i+1)=0; %boundary conditions
    u(n,i+1)=0; %boundary conditions
    v(1,i+1)=0;
    v(n,i+1)=0;

    
    if mod(i,10)==1, %Plot every 10th iterate
            clf;
        plot(x,u(:,i+1),x,v(:,i+1),'LineWidth',2);
        ylabel('$u(x,t)$','Interpreter','latex','Fontsize',18);
        xlabel('$x$','Interpreter','latex','Fontsize',18);
        title(strcat('$t=$',num2str(t(i))),'Interpreter','latex','Fontsize',18 );
        axis([0,1,0,1]);
        drawnow;
    end
    
end
