function pdex9
m = 0;
x = linspace(0,1);
t  =  linspace(0,10000000000);
sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t);
u1 = sol(:,:,1);
u2 = sol(:,:,2);
u3 = sol(:,:,3);
u4 = sol(:,:,4);
%––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––
figure
plot(x,u1(end,:))
title(u1(x,t))
xlabel(Distance x)
ylabel(u1(x,1))
%––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––
figure
plot(x,u2(end,:))
title(u2(x,t))
xlabel(Distance x)
ylabel(u2(x,2))
% –––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––
figure
plot(x,u3(end,:))
title(u3(x,t))
xlabel(Distance x)
ylabel(u3(x,3))
%––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––
figure
plot(x,u4(end,:))
title(u4(x,t))
xlabel(Distance x)
ylabel(u4(x,4))
%––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––
function [c,f,s] = pdex4pde(x,t,u,DuDx)
c = [0; 0; 0; 0];
f = [1; 1; 1; 1].* DuDx;
r1  =  0.1;r2  =  0.1; r3  =  0.1; r4  =  0.1; r5  =  0.1; a  =  0.1; b  =  0.1;
F1  = - r1  *  u(1) - (r2  *  u(1)*  u(2))/(a  *  u(1)+u(2)*  b);
F2  = - (r3  *  u(1)*  u(2))/(a  *  u(1)+u(2)*  b);
F3  = ((r4/2)*  u(1));
F4  = (r5  *  u(1)*  u(2))/(a  *  u(1)+u(2)*  b);
s  =[F1; F2; F3; F4];
% –––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––
function u0 = pdex4ic(x)                %create a initial conditions
u0 = [1; 1; 0; 0];
% –––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––
function [pl,ql,pr,qr]=pdex4bc(xl,ul,xr,ur,t)      %create a boundary conditions
pl = [0; 0; ul(3); 0];
ql = [1; 1; 0; 1];
pr = [ur(1)-1; ur(2)-1; ur(3)-0; ur(4)-0];
qr = [0; 0; 0; 0];
Algorithm 1: Matlab program to find the numerical solution of (5)–(8).