how to implement pdepe for following system of equations?

I need to solve following PDE system of equations with the outputs are sigma, rho, chi, theta. All of them are functions of (eps,y):

dsigma/deps= mu*(1-q/Q)+M d/dy(dsigma/dy);

drho/deps=alpha*(1-rho/exp(-1/chi));

dchi/deps=beta*(1-chi/chi0);

dtheta/deps=K*sigma*q/Q+(K1/Q)*d/dy(dtheta/dy)-(K2/Q)*(theta-theta0)

with following conditions:

rho(0,y)=rhoInitial; chi(0,y)=chiInitial; sigma(0,y)=sigmaInitial; theta(0,y)=theta0

sigma(eps,1)=sigma(eps,-1); theta(eps,1)=theta(eps,-1)

Does anybody have some idea to implement these equations as well as boundary conditions on Matlab? Thank you very much.

Best regards, Minh