Trying to solve a system of 5 PDE's using pdepe but getting errors?

I can't get my code to run! I get these massages:

 Attempted to access u(2); index out of bounds because numel(u)=1.

Error in chagastw>chagastwpde (line 53)
s = [(50*exp(-(1/100)*(u(1)+u(2)+u(3))))*(u(1)+qa*(1-pa)*u(2)+qc*(1-pc)*u(3))-mh*u(1)-((bvh*u(1)*u(5))/(u(1)+u(2)+u(3)))+na*u(2);...

 Error in pdepe (line 246)
[c,f,s] = feval(pde,xi(1),t(1),U,Ux,varargin{:});

 Error in chagastw (line 21)
sol = pdepe(m,@chagastwpde,@chagastwic,@chagastwbc,x,t,mh,mv,dha,dhc, qa, pa, qc,pc,na, zeta,bvh, bhav, bhcv,Dr);

My codes are

 function chagastw    
  mh=.7;
     mv=1.2; 
     dha=0.2; 
     dhc=0.5;
     qa=0.1;
     pa=0.1; 
     qc=0.3;
     pc=0.1;
     na=0.7;  
     zeta=.4; 
     bvh= 3;
     bhav=.2; 
     bhcv=2.6; 
     Dr=0.4;

 m = 0;
 x = linspace(-10,10,200);
 t = linspace(0,20,100);

 sol = pdepe(m,@chagastwpde,@chagastwic,@chagastwbc,x,t,mh,mv,dha,dhc, qa, pa, qc,pc,na, zeta,bvh, bhav, bhcv,Dr);

 u1 = sol(:,:,1);
 u2 = sol(:,:,2);
 u3 = sol(:,:,3);
 u4 = sol(:,:,4);
 u5 = sol(:,:,5);

 figure
 surf(x,t,u2) 
 title('Numerical solution of I_{ha}')
 xlabel('Distance x')
 ylabel('Time t')

 function [c,f,s] =  chagastwpde(x,t,u,DuDx,mh,mv,dha,dhc, qa, pa, qc,pc,na, zeta,bvh, bhav, bhcv,Dr)
 c = [1;1;1;1;1];
 f =[0 ;0 ;0 ;0.4 ; 0.4].*DuDx;
 s = [(50*exp(-(1/100)*(u(1)+u(2)+u(3))))*(u(1)+qa*(1-pa)*u(2)+qc*(1-pc)*u(3))-mh*u(1)- ((bvh*u(1)*u(5))/(u(1)+u(2)+u(3)))+na*u(2);...
    (50*exp(-(1/100)*(u(1)+u(2)+u(3))))*(qa*pa*u(2)+qc*pc*u(3))+((bvh*u(1)*u(5))/(u(1)+u(2)+u(3)))-(mh+dha+na+zeta)*u(2);...
    zeta*u(2)-(mh+dhc)*u(3);...
    100*(u(4)+u(5)).*exp(-(1/200)*(u(4)+u(5)))-((bhav*u(4)*u(2))/(u(1)+u(2)+u(3)))-((bhcv*u(4)*u(3))/(u(1)+u(2)+u(3)))-mv*u(4);...
    ((bhav*u(4)*u(2))/(u(1)+u(2)+u(3)))+((bhcv*u(4)*u(3))/(u(1)+u(2)+u(3)))-mv*u(5)];
% --------------------------------------------------------------
function u0= chagastwic(x,mh,mv,dha,dhc, qa, pa, qc,pc,na, zeta,bvh, bhav,bhcv,Dr);
u0 = heaviside(x+2)-heaviside(x-2);
% -------------------------------------------------------------- 
function [pl,ql,pr,qr] = chagastwbc(xl,ul,xr,ur,t,mh,mv,dha,dhc, qa, pa, qc,pc,na, zeta,bvh, bhav, bhcv,Dr)
pl = [0;0;0;0;0];
ql = [1;1;1;1;1];
pr = [0;0;0;0;0];
qr = [1;1;1;1;1];

[Bill Greene]

Possibly you want

u0 = (heaviside(x+2)-heaviside(x-2))*ones(5,1);

[AHUOD ALSHERI]

Thank you for notifying me on that.

But still the codes won't run. Now I get this massage:

Warning: Failure at t=0.000000e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (7.905050e-323) at time t. > In ode15s at 668 In pdepe at 289 In chagastw at 21 Warning: Time integration has failed. Solution is available at requested time points up to t=0.000000e+00. > In pdepe at 303 In chagastw at 21 Error using surf (line 78) Z must be a matrix, not a scalar or vector

Error in chagastw (line 31) surf(x,t,u2)

[Bill Greene]

I suggest you post the equations you are trying to solve in mathematical form.

[Torsten]

@AHOUD ALSHERI:

You divide by zero when calculating "s" in "chagastwpde".

[AHUOD ALSHERI]

This is the equations I am trying to plot