把帖子发到悬赏板块吧楼主，重赏之下必有勇夫

function ff=djf(p1,re)

pe=10;pw=4;la=0.001;lc=0.2;a=(lc-la)/lc;b=1/la;rw=0.1;

  g=@(r,p1)((b*(p1-pw)/log(re/rw)-a*r+r+sqrt((a*r-r-b*(p1-pw)/log(re/rw)).^2+4*a*b*r*(p1-pw)/log(re/rw)))./(2*b*r));
   ff=pe-pw-quadl(@(r)g(r,p1),0.1,re);
end

   
   

%主程序调用文件
% [x,feval]=fsolve(@djf,6,optimset('Display','off','TolX',1e-8,'TolFun',1e-8))
%不同生产时间，产量不同，用稳态逼近方法所对应的供给半径不同。选择不同半径，可以求得折算压力和修正系数
%以下是批处理主程序代码
% q=[1:100];
% for i=1:100
% re(i)=0.1*exp(0.0864*2*pi*10*6/q(i));
% end
% re=re';
% arrayfun(@(xx)fsolve(@(x)djf(x,xx),6),re)
%arrayfun要批处理的参数是xx.