Step Function as Input Boundary Condition to MATLAB PDEPE solver

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;

Answer 1

[Bill Greene]

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

[Siamack]

Thanks alot Bill, I used the smooth function and it worked. Just one more question and that of: your suffestion was based on the error message regarding the integration tolerance?

Regards Siamack

[Bill Greene]

The error message indicates the ODE solver wasn't able to take even a tiny, first step. So I assumed that was due to the sharp discontinuity in the BC at t=0.