function rak
%
%   The PDEs are
%   
%      D(u)/Dt = D^2(u)/Dx^2 + u(1-u-rv)
%      D(v)/Dt = D^2(v)/Dx^2 - buv
% 
%
%   In the form expected by PDEPE, the equations are
%
%   |1|         |u|      | D(u)/Dx |    | u(1-u-rv)|
%   | | .*  D_  | | = D_ |         | +  |          |
%   |1|     Dt  |v|   Dx | D(v)/Dx |    | -buv     |
%
%   ---         ---   ----------------  ---------------
%    c           u      f(x,t,u,Du/Dx)   s(x,t,u,Du/Dx)
%
%   The initial condition is u(x,0) = 0 and v(x,0) = 0 for 0 <= x <= 1.
%   The left boundary condition is D(u)/Dx = 0, v(0,t) = 1.  The 
%   condition on the partial derivative of u1 has to be written in terms 
%   of the flux.  In the form expected by PDEPE, the left bc is
%
%      |0   |       |1|     | D(u)/Dx |   |0|
%      |    |   +   | | .*  |         | = | |
%      |v -1|       |0|     | D(v)/Dx |   |0|
%
%      ---        ---     ------------------   ---
%    p(0,t,u)    q(0,t)     f(0,t,u,Du/Dx)      0
%
%   The right boundary condition is u(1,t) = 1, D(v)/Dx = 0:
%
%      |u - 1|       |0|     | D(u)/Dx |   |0|
%      |     |   +   | | .*  |         | = | |
%      |  0  |       |1|     | D(v)/Dx |   |0|
%
%      -------       -----   ----------------   ---
%      p(1,t,u)      q(1,t)  f(1,t,u,Du/Dx)      0
%
%   
m = 0;
x = linspace(0,100,200);
t = linspace(0,200,200);

sol = pdepe(m,@bz_pr,@bz_ic,@bz_bc,x,t);
u = sol(:,:,1);
v = sol(:,:,2);
%sol = pdepe(m,@bz_pr2,@bz_ic,@bz_bc,x,t);
%u2 = sol(:,:,1);
%v2 = sol(:,:,2);
% 
% figure(1);
% surf(x,t,u);
% title('u(x,t)');
% xlabel('Distance x');
% ylabel('Time t');
% 
% figure(2);
% surf(x,t,v);
% title('v(x,t)');
% xlabel('Distance x');
% ylabel('Time t');

figure(3)
for i=1:length(t)
    %subplot(2, 1, 1)
    plot(x,u(i,:),'g',x,v(i,:),'r',x,1-u(i,:),'b');
    title(t(i))
    legend('v','u')
    %ylim([0 1.2])
    %subplot(2, 1, 2)
    %plot(x,u2(i,:),'g',x,v2(i,:),'r');
    %title(t(i))
    %legend('u','v')
    %ylim([0 1.2])
    drawnow;
    
end
% --------------------------------------------------------------------------

function [c,f,s] = bz_pr(x,t,u,DuDx)

gamma=50;
thetha=1;
betha=5;
c = [1; 1];                                  
f = [0; 1] .* DuDx;                   
s1= gamma*(1-u(1)-u(1)*u(2));
s2= u(2)*(1-thetha*u(2)-betha*u(1));
s = [s1; s2]; %[0 ;0];%                                

% --------------------------------------------------------------------------
function [c,f,s] = bz_pr2(x,t,u,DuDx)

r=50;
c = [1; 1];                                  
f = [1; 1] .* DuDx;                   
s1= u(1)*(1-u(1)-r*u(2));
s2= -u(1)*u(2);
s = [0 ;0];%                                

% --------------------------------------------------------------------------
function u0 = bz_ic(x)
u0 = [1 ; exp(-(x-50)^2)/2];                                 

% --------------------------------------------------------------------------

function [pl,ql,pr,qr] = bz_bc(xl,ul,xr,ur,t)
pl = [0; 0];                               
ql = [1; 1];                                  
pr = [0; 0];                            
qr = [1; 1];                                  

