BioHeat Equation solution in MATLAB using pdepe

[m.r.fouladi]

We have this Equation as bioheat equation:
∂T/∂t = α ∇²T + 1/ρc[S+S_p+S_m]
and also this:
S_p=m_bc_b(T_ab-T)
that all α,ρ,c,S,S_m,m_b,c_b,T_ab are constants, now I want to solve this equation in conditions below with pdepe in MATLAB:
There is a Tumor as a sphere with radius 1 cm exactly in center of a Normal Tissue with radius of 5 cm, an electrode at t=0 gives an Energy to the center of Tumor for 400 seconds that T will vary and we want to have this variation by solving this equation, also below assumptions is considered:
We assume that Tumor has the radius of 1 cm and that is perfused at the same level as healthy tissue. The tumor is located at the center of a sphere of healthy tissue that has a radius of 5 cm and all of the tissue has a metabolic energy release rate of 145 W/m³.
The power is delivered locally to the tumor by positioning the electrode in the center of the tumor. the initial temperature of tumor and surrounding tissue is uniform at 37 °C. in addition, It is assumed that because of symmetry there is zero flux at the center of the tumor and that at the outer surface of healthy tissue the flux is also zero. the following numerical values for the parameters are used:
m_b = 0.18
c_b = 3300
S_m = 145
T_ab = 37 °C
α = 10∧(-7)
c=3800
S= 4 * 10∧(5)
ρ= 850.I coded it in MATLAB such below, but my code isn't right as my knowledge, please correct it my masters:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function pdex1

m = 2;
x = 0.001:0.001:05;
t = 0:1:400;

sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
% Extract the first solution component as u.
u = sol(:,:,1);

% A surface plot is often a good way to study a solution.
surf(x,t,u)
title('Numerical solution computed with 20 mesh points.')
xlabel('Distance x')
ylabel('Time t')

% A solution profile can also be illuminating.
figure
plot(x,u(end,:))
title('Solution at t = 2')
xlabel('Distance x')
ylabel('u(x,400)')
% --------------------------------------------------------------
function [c,f,s] = pdex1pde(x,t,u,DuDx)
c = 1e7;
f = DuDx;
if (x<0.01)
s = (18.089*10^(-2)-(18*10^(-5)*u))*(1e7);
else
s=594*(310-u)/(850*3800)
end
% --------------------------------------------------------------
function u0 = pdex1ic(x)
u0 = 310;
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
pl = ul-exp((-18)*10^(-5)*t)-1004.94;
ql = 0;
pr = ul-exp((-18)*10^(-5)*t)-1004.94;
qr = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

why It's failed?

 

[DrClaude]

why It's failed?

Why do you say it failed?

Also, I have added [code] [/code] tags around your code to make it more readable. Please do that every time you post code.

 

[m.r.fouladi]

Why do you say it failed?

Also, I have added [code] [/code] tags around your code to make it more readable. Please do that every time you post code.

thanks for your attention.
because It doesn't agree with my physics knowledge of this problem