I'm trying to solve the following system of differential equations numerically. What are the available finite difference approaches and matlab solvers to solve such a system? Other approaches to solve the problem are welcome as well. The first two equations are quasi-steady and the rest of the equations have time derivatives in them. I'm trying using matlab's pdepe solver, but since the boundary conditions have to be specified even for the ODEs which have time-derivatives, I didn't think it as a viable option.
$$0= \rho \frac{\partial}{\partial z}\left(D \frac{\partial c_1}{\partial z} \right) - \it{source_1}\, $$
$$ 0= \rho \frac{\partial}{\partial z}\left(D \frac{\partial c_2}{\partial z} \right) - \it{source_2}\,$$
$$\frac{\partial c_3}{\partial t} = source_3$$
$$\frac{\partial c_4}{\partial t} = source_4$$
$$D = A_1 c_1 + A_2 c_2 + \dots + A_n c_n \,, \rho = 1 (say)\,.$$
The source terms for any particular equation (or concentration) contain dependence on the other concentrations.