function pdex4
global u1
global u2
global x
      
m = 0;
x = 0:1:900; % size of pond under consideration
tend =2000;  % time till computation is performed
t = 0:1:tend;

sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t);
u1 = sol(:,:,1); %Prey
u2 = sol(:,:,2); %Predator

% --------------------------------------------------------------
function [c,f,s] = pdex4pde(x,t,u,DuDx)
% subfunction 1, defines the partial differential equation
% parameters chosen as given in the paper
h = 0.3;
k = 2.0;
r = 0.4;
m = r*k;
c = [1; 1];                                  
f = [1; 1] .* DuDx;                   
F1 = u(1)*(1-u(1)) - (u(1)*u(2))/(u(1) + h); % F1 = u(1-u) - uv/(u+h)
F2 = (k*u(1)*u(2))/(u(1)+h) - m*u(2);        % F2 = kuv/(u+h) - mv
s = [F1; F2];                                                          
% --------------------------------------------------------------
function u0 = pdex4ic(x)
% subfunction 2, defines the initial conditions
% conditions similar to those given in the paper
r = 0.4;
h = 0.3;
% disp(x)
epsilon = 10^-2;  %10^-6 is the -4.5
delta = -1.5*10^-3;
ustar =r*h/(1-r);
vstar = (1-ustar)*(h+ustar);
% u0 = [ustar; (vstar + delta+ epsilon.*x) ];   for non uniform initial
% conditons
A = 5;
S = 400;
xo=600;
if x >=50 && x<=60
    a=0.3;
u0 = [ustar; (vstar + delta+ epsilon.*x + a)];
else
u0 = [ustar; (vstar + delta+ epsilon.*x)]     ;                            
end
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex4bc(xl,ul,xr,ur,t)
% subfunction 3, defines the boundary conditions
% Neuman zero flux boundary conditions used
pl = [0; 0];                               
ql = [1; 1];                                  
pr = [0; 0];                            
qr = [1; 1];           