pdepe with discontinuos domain

Hello All,

I am trying to implement a set of pdes (1-D, time dependent) having parameters with different values in different domains (for example, 0<x<17, 17<x<34). The interfaces of domain are total 6 in number. The code is working fine (without errors).

Problem: The obtained solution is not smooth at the given interface, although a lot of points of domain are defined near the interface.

The code is as follows. Any kind of help is appreciated.

Thanks and regards,

-Meenal

function abc

clear all; clc;
global Temperature L tend cp k rho wb qm cb  Ta  rho_b qs 
L=0.093; 
tend=3; 
Ta=37; 
cp=[3800; 3700; 3800; 3700; 3800; 2300; 4300]; 
k=[0.5; 0.5; 0.5; 0.5; 0.5; 1.16; 0.34]; 
rho=[1052; 1052; 1052; 1052; 1052; 1052; 1000]; 
wb=[0.00833; 0.0005; 0.00833; 0.0005; 0.00833; 0.0003; 0.00033]; 
qm=[25000; 25000; 25000; 25000; 25000; 3700; 3600]; 
cb=3800; 
rho_b=1052;
qs=0;
m=0;

x = [linspace(0,17,40) linspace(17.01,34,40) linspace(34.01,51,40) linspace(51.01,68,40) linspace(68.01,85,40) linspace(85.01,89,10) linspace(89.01,93,10)];
t=linspace(0,tend,41);

sol = pdepe(m,@pdex2pde,@pdex2ic,@pdex2bc,x,t);

Temperature = sol(:,:,1);
figure
surf(x,t,Temperature)
title('Numerical solution with nonuniform mesh')
xlabel('Distance x')
ylabel('Time t')
zlabel('Solution Temperature')
figure
plot(x,Temperature(end,:))
xlabel('Distance x')
ylabel('Solution Temperature')
title('Solution profiles at several times')



function [c,f,s] = pdex2pde(x,t,T,DTDx)
global cp k rho wb qm cb  Ta  rho_b qs 
if x <= 17
    c = rho(1).*cp(1);
    f = k(1).*(DTDx);
    s = wb(1).*cb*rho_b.*(Ta-T)+qm(1)+qs;
elseif (x>17 && x<=34)
   c = rho(2).*cp(2);
    f = k(2).*(DTDx);
    s = wb(2).*cb*rho_b.*(Ta-T)+qm(2)+qs;
elseif (x>34 && x<=51)
    c = rho(3).*cp(3);
    f = k(3).*(DTDx);
    s = wb(3).*cb*rho_b.*(Ta-T)+qm(3)+qs;
elseif (x>51 && x<=68)
    c = rho(4).*cp(4);
    f = k(4).*(DTDx);
    s = wb(4).*cb*rho_b.*(Ta-T)+qm(4)+qs;
elseif (x>68 && x<=85)
    c = rho(5).*cp(5);
    f = k(5).*(DTDx);
    s = wb(5).*cb*rho_b.*(Ta-T)+qm(5)+qs;
elseif (x>85 && x<=89)
    c = rho(6).*cp(6);
    f = k(6).*(DTDx);
    s = wb(6).*cb*rho_b.*(Ta-T)+qm(6)+qs;
elseif (x>89 && x<=93)
    c = rho(7).*cp(7);
    f = k(7).*(DTDx);
    s = wb(7).*cb*rho_b.*(Ta-T)+qm(7)+qs;
end

function T0 = pdex2ic(x)
T0 =309;

function [pl,ql,pr,qr] = pdex2bc(xl,Tl,xr,Tr,t)
pl = 0; 
ql = 1; 
pr = Tr - 312;
qr = 0;