# How to solve these 2-D coupled differential equations in Matlab?

I have three partial differential equations (PDEs) and an analytical solution for a variable as shown. Using these equations I want to solve for \phi(x,y,t), p(x,y,t), C_{a}(x,y,t) and C_{b}(x,y,t) i.e. in terms of space and time.

I know there is a function pdepe( ) in Matlab to solve initial-boundary value problems for parabolic-elliptic PDEs in 1-D. I would like to know how this function or some other in Matlab can be used to solve the problem described below which is 2-D and coupled.

PROBLEM:

Following two equations represent PDEs for two species a and b, respectively:

Where D_{h} and q are given as:

Here, R_{a}=R_{b}=R, where R is given as:

Finally, the last equation is given as:

INITIAL AND BOUNDARY CONDITIONS:

Total domain size is 10 cm x 5 cm and the width of the y-shaped subdomain is 0.5 cm. This subdomain has an initial \phi of 0.50 while in the surrounding matrix \phi= 0.26. Constant p of 1 Pa and 0 Pa are maintained at boundaries (1) and (2) respectively, corresponding to a gradient of approximately 10^-3 m m^-1. The p on boundaries (3) and (4) are determined by linear gradients between boundaries (1) and (2). Constant C of C_{a} = 2 mol m^-3 and C_{b} = 0.2302 mol m^-3 are maintained at boundary (3), while the concentrations at boundary (4) are set at C_{a} = 1 mol m^-3 and C_{b} = 0.4603 mol m^-3. The concentrations at boundary (1) are determined by constant gradients between boundaries (3) and (4), while an advective flux boundary condition $$(\frac{\partial C}{\partial x} = 0)$$ is set at the outlet at (2).

-