MATLAB: Step Function as Input Boundary Condition to MATLAB PDEPE solver

1d heat conductionmatlab pdepe toolbox

Dear Matlab Users I have the following code to find the transient solution to the 1D heat conduction equation. The code below intends to use a rectangular heat pulse as the initial temperature but i receive an error ( Warning: Failure at t=0.000000e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (7.905050e-323) at time t). Can anyone tell me what is issue with the boundary condition as I set it and if there is a better way to deal with the matter.
function out = parabolic_silicon(~)
global rho cp k T_i T_o
L   = 10;            % [cm]
k   = 1.3;          % [W/cmK]
rho = 2.33;         % [g/cm^3]
cp  = 0.7;          % [J/gK]
T_i = 295;          % [K]
T_o = T_i + 5.2e-3;  % [K] 
t_end= 1;     % [s]
m = 0;
x = linspace(0,L,500);
t = linspace(0,t_end,500);
sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
% Extract the first solution component as u.
Temperature = sol(:,:,1);
L = strcat('t = ',num2str(t(:)),' Seconds');
plot(x,sol);
% legend(L,'location','northeast');
% % --------------------------------------------------------------

% OUTPUT
out = {x,t,sol};
% % --------------------------------------------------------------
function [c,f,s] = pdex1pde(x,t,u,DuDx)
global rho cp k
c = rho*cp;
f = k*DuDx;
s = 0;
% --------------------------------------------------------------

function u0 = pdex1ic(x)
u0 = 295;
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
global T_i T_o
pl = ul - (T_i + (T_o*(heaviside(t) - heaviside(t - 0.5))));
ql = 0;
pr = ur - T_i;
qr = 0;

Best Answer

I suggest you smooth the sharp corners in your boundary condition profile using a Heaviside function similar to the code below. Note, you can change the sharpness of the corners by varying the r variable.
function k=heavisideSmooth(t)
  k1 = 0; k2 = 1; c = 0;
  r = 1e3;
  ecr = exp(c*r);
  erx = exp(r*t);
  k = (k1*ecr + k2*erx)./(ecr+erx);
end
Related Question