function Competition1D
%
% Evolution of the competition system under linear-quadratic interaction:
%
% y_t - K_1 y_{xx} = Ay - Byz
% z_t - K_1 z_{xx} = Cz - Dyz
%
% together with Dirichlet and initial conditions
%

format short;
% format long;
close all;

global K1 K2 A B C D IPLOT NX1 T Y0 NT1 XNUMER HX NX

%
% THE DATA OF THE PROBLEM
%

K1=.1; K2=.1;                  % DIFFUSION COEFFICIENTS
A=.5; B=.08; C=1.1; D=1.;        % BIRTH/DEATH AND INTERACTION COEFFICIENTS

T=10.;                          % FINAL TIME
Lmax=1.;                        % HALF-LENGTH

y1max=.05; sigma1=.05;         % DATA FOR INITIAL PREY POPULATION
y2max=.1; sigma2=.15;         % DATA FOR INITIAL PREY POPULATION

NX=50;                         % FOR SPATIAL APPROXIMATION
% NX=20;                         % FOR SPATIAL APPROXIMATION
NT=100;                         % NUMBER OF TIME STEPS
% NT=20;                         % NUMBER OF TIME STEPS

IPLOT=2;                       % FIGURE REPRESENTATION LEVEL
% IDISP=1;                       % PRINT LEVEL

%
% PREPARING COMPUTATIONS
%

HX=2.*Lmax/NX;                  % SPACE STEP
NX1=NX+1; 
XNUMER=linspace(-Lmax,Lmax,NX1);% THE x_j

% HT=T/Ntt;

coef1=.5/sigma1; coef2=.5/sigma2;                
Y0=[y1max*exp(-coef1*XNUMER.*XNUMER);y2max*exp(-coef2*XNUMER.*XNUMER)];

NT1=NT+1;

%

tnumer=linspace(0.,T,NT1);
sol=pdepe(0,@PDEfun,@ICfun,@BCfun,XNUMER,tnumer);
    
if (IPLOT > 0)
    figure(1);
    surf(XNUMER,tnumer,sol(:,:,1));
    xlabel('x: SPACE');
    ylabel('t: TIME');
    zlabel('y: FIRST SPECIES');

    figure(2);
    surf(XNUMER,tnumer,sol(:,:,2));
    xlabel('x: SPACE');
    ylabel('t: TIME');
    zlabel('z: SECOND SPECIES');
end


% plot(x, y, 'k', x, y, 'r*')
% legend('Interpolación de Lagrange',...
%        'Interpolación lineal a trozos',...
%        'Location', 'Best')
% 
% format long;

shg;

end

function [c,f,s]=PDEfun(~,~,u,Du)

% [C,F,S]=PDEFUN(X,T,U,DUDX)

global K1 K2 A B C D IPLOT NX1 T Y0 NT1 XNUMER HX NX

c=[1.;1.];
f=[K1*Du(1);K2*Du(2)];
s=[A*u(1) - B*u(1)*u(2);C*u(2) - D*u(1)*u(2)];

end

function [u]=ICfun(x)

%

global K1 K2 A B C D IPLOT NX1 T Y0 NT1 XNUMER HX NX

u=interp1(XNUMER,Y0',x);

end

function [pl,ql,pr,qr]=BCfun(~,ul,~,ur,~)

% [PL,QL,PR,QR]=BCFUN(XL,UL,XR,UR,T)

pl=ul; ql=[0.;0.];
pr=ur; qr=[0.;0.];

end
