I have the following system of equations which I'm trying to solve using Matlab's pdepe solver.The 1-D spherical heat diffusion equation with heat generation (source term):$$\rho \, C_p\frac{\partial T}{\partial t} = \frac{1}{r^2} \frac{\partial}{\partial r} \left(r^2 \, k \frac{\partial T}{\partial r} \right) + \Delta H \frac{\partial \rho_w}{\partial t} \\$$and boundary conditions:$$k \frac{\partial T}{\partial r} = h\,(T_\infty - T_s) \; \; \text{at r = surface}\\\frac{\partial T}{\partial r} = 0 \; \; \text{at r = 0}$$also note that $\Delta H$ is a constant, i.e. $\Delta H = -255,000$. The initial value for $\rho_w = \rho_{initial} = 700$ and the ambient temperature $T_\infty = 773$ (see Matlab code below for more details).The kinetic reactions of the system:$$\frac{\partial \rho_w}{\partial t} = -(K_1+K_2+K_3) \rho_w \\\frac{\partial \rho_g}{\partial t} = K_1 \rho_w \\\frac{\partial \rho_T}{\partial t} = K_2 \rho_w \\\frac{\partial \rho_c}{\partial t} = K_3 \rho_w \\\frac{\partial \rho_g}{\partial t} = K_4 \rho_T \\\frac{\partial \rho_c}{\partial t} = K_5 \rho_T \\\text{where}\; \; K = A\,e^{\frac{-E}{RT}}$$The thermal parameters $\rho, C_p, k$ vary according to the following:$$\rho = \rho_w + \rho_c \\Y_w = \frac{\rho_w}{\rho_{initial}} \\C_p = Y_w C_{pw} + (1-Y_w) C_{pc} \\k = Y_w k_w + (1-Y_w) k_c $$where $C_{pw} = 1500$, $C_{pc} = 1100$, $k_w = 0.105$, and $k_c = 0.071$.I can solve the single PDE for heat diffusion (no source term) using the pdepe function in Matlab as follows:function PDEexampleSphered350 = 0.035e-2; % diameterr350 = d350/2;   % radiusTi = 300;Tinf = 773;tmax = 0.8;m = 2;x = linspace(0,r350,20);t = 0:0.01:tmax;sol = pdepe(m,@pdefunc,@icfunc,@bcfunc,x,t);u = sol;% surface plotfigure(3)surf(x,t,u) xlabel('Distance (m)')ylabel('Time (s)')% temperature profilefigure(4)plot(t,u(:,1),'b',t,u(:,end),'r')hold onplot([0 tmax],[Tinf Tinf],':k')hold offaxis([0 tmax Ti-20 Tinf+20])xlabel('Time (s)')ylabel('Temperature (K)')% --------------------------------------------------------------function [c,f,s] = pdefunc(x,t,u,dudx)rho = 700;  % densitycp = 1500;  % heat capacityk = 0.105;  % thermal conductivityc = rho*cp;f = k.*dudx;s = 0;% --------------------------------------------------------------function u0 = icfunc(x)Ti = 300;   % initial temperatureu0 = Ti;% --------------------------------------------------------------function [pl,ql,pr,qr] = bcfunc(xl,ul,xr,ur,t)h = 375;    % heat transfer coefficientTinf = 773; % ambient temperaturepl = 0;ql = 0;pr = -h*(Tinf-ur);qr = 1;How can I include the kinetic reactions and heat generation term into the Matlab function?Note - for more information on the pdepe function in Matlab, please see the following documentation: http://www.mathworks.com/help/matlab/ref/pdepe.html
How to add reaction and source terms to a diffusion PDE solver written with MATLAB's pdepe function?
pde;matlab;finite difference;parabolic pde;implicit methods