%% MATLAB CODE FOR SOLVING THE PDE MODEL IN  
%% THE PROJECT ON INVASION OF DESERTS. 
%% WRITTEN BY JONATHAN SHERRATT.
clear all

%% Set equation parameters
aparameter=1.5;
bparameter=0.45;
nuparameter=0; % LEAVE THIS AS ZERO UNTIL THE FINAL PART OF THE LAB

%% Set spatial and temporal resolution parameters
xmax=150;
npts=301;
dx=xmax/(npts-1);
tmax=60;
dt=0.01;
theta=0.5;
dtodx2=dt/(dx*dx);
dtnuodx=nuparameter*dt/dx;

%%set output parameters
nplot=8;
niters=(tmax/dt)/nplot;

uout=zeros(nplot,npts);
wout=zeros(nplot,npts);
tout=zeros(nplot,1);
timerec=0;

%% Initialise space coordinates
x=[0:dx:xmax];

%% Initialise u and w
us = (aparameter+sqrt(aparameter*aparameter-4*(bparameter*bparameter)))/(2*bparameter);
ws = bparameter/us;
u=zeros(npts,1);
w=zeros(npts,1);
for jj=1:npts
  if jj>(npts/2) 
    u(jj) = 0;
    w(jj) = aparameter;
  else
    u(jj) = us;
    w(jj) = ws;
  end
end
u0=u;
w0=w;

%% Set up a sparse matrix for diffusion that is always used
A=sparse(npts*2,npts*2);

for jj=2:(npts-1)
    A(jj,jj)=1+2*dtodx2*theta;
    A(jj+npts,jj+npts)=1;
    A(jj,jj+1)=-dtodx2*theta;
    A(jj,jj-1)=-dtodx2*theta;
end

%% Set boundary conditions
A(1,1)=1;
A(npts,npts)=1;
A(npts+1,npts+1)=1;
A(2*npts,2*npts)=1;

%% Misc other setups
cc=[2:(npts-1)];
udiff=zeros(npts,1);
wdiff=zeros(npts,1);
ukinetics=zeros(npts,1);
wkinetics=zeros(npts,1);
b=zeros(npts*2,1);

%% Main solution loop

for iplot=1:nplot %For each time up to when you want data output
    for tt=1:niters %for each iteration within that

       timerec=timerec+dt;
        
       for jj=2:(npts-1)
         ukinetics(jj,1)=(u(jj,1)*u(jj,1)*w(jj,1)-bparameter*u(jj,1))*dt;
         wkinetics(jj,1)=(aparameter-w(jj)-u(jj)*u(jj)*w(jj))*dt;
       end
       
       udiff(cc,1)=(u(cc+1)-2*u(cc)+u(cc-1))*dtodx2*(1-theta);
       if nuparameter>0.0
         wdiff(cc,1)=(w(cc+1)-w(cc))*dtnuodx;
       else
         wdiff(cc,1)=(w(cc)-w(cc-1))*dtnuodx;
       end
       
       b(cc,1)=u(cc)+udiff(cc)+ukinetics(cc);
       b(cc+npts,1)=w(cc)+wdiff(cc)+wkinetics(cc);
       b(1,1)=us;
       b(npts,1)=0;
       b(npts+1,1)=ws;
       b(2*npts,1)=aparameter;

       bb=A\b;
       
       u=bb(1:npts,1);
       w=bb(npts+1:2*npts,1);
    end
    uout(iplot,:)=u;
    wout(iplot,:)=w;
    tout(iplot,1)=timerec;
end

%% End of main solution loop

%% Plot the results

figure;
subplot(2,1,1)
plot(x,u0,x,uout(1,:),x,uout(2,:),x,uout(3,:),x,uout(4,:),x,uout(5,:),x,uout(6,:),x,uout(7,:),x,uout(8,:))
ylabel('Vegetation biomass, u(x)');
xlabel('Space, x');
title(strcat('A=',num2str(aparameter),'  B=',num2str(bparameter),'  nu=',num2str(nuparameter),'  tmax=',num2str(tmax)));
subplot(2,1,2)
plot(x,w0,x,wout(1,:),x,wout(2,:),x,wout(3,:),x,wout(4,:),x,wout(5,:),x,wout(6,:),x,wout(7,:),x,wout(8,:))
ylabel('Water density, w(x)');
xlabel('Space, x');

%% End of program
