function Conidiation1D
%
% Evolution of the 1D conidiation system:
%
% e_t - D e_{xx} = F(e) + u k
%
% e = (e_1,...e_6)^t, D = diag(d_i),
% F(e) = (0,                   a_{12}e_1, ...
%         a_{13}e_1+a_{23}e_2, 0, ...
%        -b_{15}e_1e_5,        a_{36}e_3+a_{46}e_4-b_{56}e_5e_6)
% u = u(t) (scalar), k =(1,1,1,1,0,0)^t
%
% together with Neumann and initial conditions
%

format short;
% format long;
close all;

global D A12 A13 A23 A36 A46 B15 B56 UL0 IPLOT NX1 T TS2 TS4 ...
    Y0 NT1 XNUMER HX NX

%
% THE DATA OF THE PROBLEM
%

D = [.1,.1,.1,.1,.05,2.]';

A12 = 1.;
A13 = 1.;
A23 = 1.;
A36 = .1;
A46 = .1;

B15 = .01;
B56 = 50.;

UL0 = 1.;

T=10.;                          % FINAL TIME

TS2=.5*T;
TS4=.25*T;

Lmax=1.;                        % HALF-LENGTH

e5max=10.; sigma5=.05;

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;

coef5=.5/sigma5;

Y0 = zeros(6,NX1);
Y0(5,:) = e5max*exp(-coef5*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 EZ');

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

    figure(3);
    surf(XNUMER,tnumer,sol(:,:,3));
    xlabel('x: SPACE');
    ylabel('t: TIME');
    zlabel('z: THIRD EZ');

    figure(4);
    surf(XNUMER,tnumer,sol(:,:,4));
    xlabel('x: SPACE');
    ylabel('t: TIME');
    zlabel('z: FOURTH EZ');

    figure(5);
    surf(XNUMER,tnumer,sol(:,:,5));
    xlabel('x: SPACE');
    ylabel('t: TIME');
    zlabel('z: STOPPING VARIABLE');

    figure(6);
    surf(XNUMER,tnumer,sol(:,:,6));
    xlabel('x: SPACE');
    ylabel('t: TIME');
    zlabel('z: FINAL PRODUCT');

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(~,t,u,Du)

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

global D A12 A13 A23 A36 A46 B15 B56 UL0 IPLOT NX1 T TS2 TS4 ...
    Y0 NT1 XNUMER HX NX

c = [1.;1.;1.;1.;1.;1.];
f = D.*Du;

ult = light(t);
s = [ult; A12*u(1)+ult;   A13*u(1)+A23*u(2)+ult; ...
     ult; -B15*u(1)*u(5); A36*u(3)+A46*u(4)-B56*u(5)*u(6)];

end

function u = ICfun(x)

%

global D A12 A13 A23 A36 A46 B15 B56 UL0 IPLOT NX1 T TS2 TS4 ...
    Y0 NT1 XNUMER HX NX

u = zeros(6,1);
u(5) = interp1(XNUMER,Y0(5,:),x);

end

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

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

global D A12 A13 A23 A36 A46 B15 B56 UL0 IPLOT NX1 T TS2 TS4 ...
    Y0 NT1 XNUMER HX NX

pl=ul; ql=zeros(size(ul));
pr=ur; qr=zeros(size(ur));

% pl=zeros(size(ul)); ql=D;
% pr=zeros(size(ur)); qr=D;

end


function ulight = light(t)

%

global D A12 A13 A23 A36 A46 B15 B56 UL0 IPLOT NX1 T TS2 TS4 ...
    Y0 NT1 XNUMER HX NX

ulight = 0.;

if (TS4 <= t <= TS2)
    
    ulight = UL0;

end

end



