 function shi1
p=0;
l=5;
x=linspace(0,l*pi,25);
t=linspace(0,50,500);

sol = pdepe(p,@shi1pde,@shi1ic,@shi1bc,x,t);
u1 = sol(:,:,1);
u2 = sol(:,:,2);

figure
surf(x,t,u1)
title('u1(x,t)')
xlabel('Distance x')
ylabel('Time t')

figure
surf(x,t,u2)
title('u2(x,t)')
xlabel('Distance x')
ylabel('Time t')
% --------------------------------------------------------------
function [c,f,s] = shi1pde(x,t,u,DuDx) %reaction-diffusion equations
alpha=6; d=0.5; m=2; l=5; %parameters
 h1=5*alpha-u(1)-4*u(1)*u(2)/(1+u(1)*u(1)); %equation 1
 h2=m*u(1)-m*u(1)*u(2)/(1+u(1)*u(1)); %equation 2
c = [1; 1];  %coefficient of time derivatives, usually just set as 1                             
f = [1; d*m] .* DuDx;  %diffusion part
s = [h1; h2];  %reaction part                                                
% --------------------------------------------------------------
function u0 = shi1ic(x); %initial condition
u0 = [6+0.6*cos(x);(6*6+1)+0.5*cos(x)];                                 
% --------------------------------------------------------------
function [pl,ql,pr,qr] = shi1bc(xl,ul,xr,ur,t)% boundary condition
pl = [0; 0];                               
ql = [1; 1];                                  
pr = [0; 0];                            
qr = [1; 1];        