如何修改pdepe.m函数中的边界条件

如何把pdepe 边界条件的函数：[pL,qL,pR,qR] = bcfun (xL, uL, xR, uR, t)改变成[pL,qL,pR,qR]=bcfun(xL,uL,xR,uR,u,t)，论文中给出的修改方法不理解
[pL,qL,pR,qR]=feval(bc,xmesh(1),y0(:,1),xmesh(nx),y0(:,n
x),y0,t(1),varargin{:});
[pL,qL,pR,qR]=feval(bc,x(1),u(:,1),x(nx),u(:,nx),u,tnow,var
argin{:});

这个意思很好理解呢。
matlab中，feval函数的作用：
使用函数的名称或其句柄以及输入参数 x1,...,xM 来计算函数的结果。
[y1,...,yN] = feval(fun,x1,...,xM)
feval表示的是，调用函数fun，然后输入参数分别是x1，x2，……xM

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

这句话的意思是，调用函数bc，然后bc函数所需的参数列表是xmesh(1), y0(:,1), xmesh(nx), y0(:,nx), y0, t(1), varargin{:}七个参数，最后一个参数varargin{:}指的是其他可能会输入的参数，其他可能会输入的参数不一定有，但是有的话也能处理。然后输出参数是四个pL,qL,pR,qR

[pL,qL,pR,qR]=feval(bc,x(1),u(:,1),x(nx),u(:,nx),u,tnow,var
argin{:});

这句话的意思是，同样调用函数bc，然后bc函数所需的参数列表是x(1), u(:,1), x(nx), u(:,nx), u, tnow, varargin{:}七个参数，最后一个参数varargin{:}指的是其他可能会输入的参数，同样的，其他可能会输入的参数不一定有，但是有的话也能处理。然后最后输出参数也是四个pL,qL,pR,qR
其实就是将输入参数由原来的xmesh(1), y0(:,1), xmesh(nx), y0(:,nx), y0, t(1), varargin{:}改成现在的x(1), u(:,1), x(nx), u(:,nx), u, tnow, varargin{:}，你可以看到每个参数位置上参数的变化。
以上是我对你说的语法的解释，有帮助还望给个采纳支持一下答主答题，问题内有疑问可以继续探讨哦

sol=pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t)
m=2;
x = linspace(0,1,50);
t = linspace(0.1,pi,50);
[pL,qL,pR,qR]=feval(bc,xmesh(1),y0(:,1),xmesh(nx),y0(:,nx),y0,t(1),varargin{:});

[pL,qL,pR,qR]=feval(bc,x(1),u(:,1),x(nx),u(:,nx),u,tnow,varargin{:});

surf(x,t,sol(:,:,1))

figure;

kdf=[d0,kf];
function [c,f,s]=pdex4pde(x,t,u,dudx) 
global kdf 
c=1/(kdf(1))/exp(1*u./(0.737*0.2079^(1/1.41))); 
f=dudx;
s=(1/(0.737*0.2079^(1/1.41)))*(dudx.^2); 
end

function u0=pdex4ic(x)
u0=0; 
end


function [pl,ql,pr,qr]=pdex4bc(xl,ul,xr,ur,u,t)


global kdf 
u=abs(u);
pl=0;
ql=0;
n=1000; 
r=0.000593/2; 
                    
h=r/(n-1);
s=zeros (1000,1);
a=0; 
for i=1:1:n 
    if i==1 
      s(i,1)=4*pi/3*(h/2)^3.*u(:,i);
      a=a+s(i,1); 
    end 
    if i==n 
      s(i,1)=(4*pi/3*r^3-4*pi/3*(r-h/2)^3).*u(:,i);
      a=a+s(i,1); 
    end 
    if i>1 & i<n    
     s(i,1)=(4*pi/3*(h*(i-1)+h/2)^3-4*pi/3*(h*(i-1)-h/2)^3).* u(:,i);
     a=a+s(i,1); 
    end 
end 
ss=(100*0.2079-30*a/(4/3*pi*(r)^3))/(100-30*a/(4/3*pi*(r)^3))-(u(:,n)/0.737)^1.41;  
pr=kdf(2)*ss;
qr=-kdf(1)*(exp(1*ul/(0.737*0.2079^(1/1.41))));
end