How to define a step function as input BC in pdepe?

I'm trying to define a stepchange of inflow concentration(cin) in my system. The following is my code, in line 27 cin should be from 10 to 0 at t=50, so I defined cin = 10.*heaviside(50-t), but there is always an error :

Error using vertcat

Dimensions of arrays being concatenated are not consistent.

Error in ADSmainsolver>slowsorpbc (line 53)

pl = [ul(1)-cin;0];

Error in pdepe (line 250)

[pL,qL,pR,qR] = feval(bc,xmesh(1),y0(:,1),xmesh(nx),y0(:,nx),t(1),varargin{:});

Error in ADSmainsolver (line 30)

c = pdepe(0,@transpdes,@slowsorpic,@slowsorpbc,x,t,options,D,v,theta,rhob,F,kd,k,c0,cin);

Here's my code, could anyone help me,plzzzzzz :

%ADSmainsolver

clear all

clc

close all

% numerical PDE solution for transport (a+d) and nonideal sorption

%--------------------------------------------------------------------------

T = 156.5; % maximum time [min]

T1 = 1.5; % minimum time [min]

L = 15; % length [m]

D = 0.4; % diffusivity [cm^2/min]

v = 1; % real fluid velocity [cm/min]

theta = 0.4; % porosity

rhob = 1.6; % porous medium bulk density [g/cm^3]

c0f = 0; % initial concentration in fluid [ug/L]

c0s = 0; % initial concentration in solid [ug/kg]

c0=[c0f;c0s];

F = 1; %fraction of ideal partition

kd = 1; %partition coefficient [kg/L]

k = 0; % sorption rate limit const. [1/min]

M = 32; % number of timesteps (>2)

N = 100; % number of nodes

%-------------------------- output parameters------------------------------

gplot = 1; % =1: breakthrough curves; =2: profiles 1

t = linspace (T1,T,M); % time discretization

x = linspace (0,L,N); % space discretization

cin = 10.*heaviside(50-t);% inflow concentration [ug/L]

%----------------------solver-------------------------------------------

options = odeset;

c = pdepe(0,@transpdes,@slowsorpic,@slowsorpbc,x,t,options,D,v,theta,rhob,F,kd,k,c0,cin);

Y=c(:,N,1);

%---------------------- graphical output ----------------------------------

switch gplot

case 1

plot (t,c(:,N,1)) % breakthrough curves

xlabel ('time'); ylabel ('concentration');

case 2

plot (x,c(:,:,2)','--') % profiles

xlabel ('space'); ylabel ('concentration');

end

% --------------------------------------------------------------------------

function [c,f,s] = transpdes(x,t,u,DuDx,D,v,theta,rhob,F,kd,k,c0,cin)

c = [1+rhob*F*kd/theta;1];

f = [D;0].*DuDx;

s = -[v;0].*DuDx - ([k*(1-F)*kd,-k]*u)*[rhob/theta;-1];

end

function u0 = slowsorpic(x,D,v,theta,rhob,F,kd,k,c0,cin)

u0 = c0;

end

% --------------------------------------------------------------------------

function [pl,ql,pr,qr] = slowsorpbc(xl,ul,xr,ur,t,D,v,theta,rhob,F,kd,k,c0,cin)

pl = [ul(1)-cin;0];

ql = [0;1];

pr = [0;0];

qr = [1;1];

end

Answer 1

[Rohit Pappu]

There’s an error in Line 27 because, cin is of size 1x32 and MATLAB cannot vertically concatanate it with a scalar.

To define a unit step function, without using Heaviside function :

unitstep = 10*ones(size(t)); 
unitstep(t>=50) = 0; 

Answer 2

[Bill Greene]

Here is the way you want to define such a BC:

if(t>=50)
  pl = [ul(1);0];
else
  pl = [ul(1)-10;0];
end

One thing you have to remember is that boundary condition function may be called for any value of time between the starting and ending times-- not just the times where you have requested output.