How do i solve multidimensional coupled PDE for Extraction (Shrinking Core Model)

Hello, im doing a research about kinetics study of lipids extraction and want to use shrinking core model for that.

so basically there's 4 mass transport equation (one for reactor, three for particle) needed to be solve to get profile of Cb (concentration in reactor porous) over time. the variables here are Cb, Cp, q, and rc.

the first attachment is 4 equation for mass transport using shrinking core model (the first one is for reactor scale, the next three is for particle scale). the second attachment is the boundary condition.

is there any recommendation on how to solve this problem? or what method do i have to use? it will be a great help if you could provide a code to solve this kind of problem