希望大家可以matlab下，我自个搞了好几天，问了同学，还是没出来。

首先，将所有的已知参数赋值，留下仅包含未知变量的积分方程；
 
其次，将积分方程用某种积分方法如辛普森方法或其它方法离散化，得到具有特定步长的离散方程组；
 
最后，解出方程组的解，或让某个变量在某个范围内变化，求出其它变量的变化结果并绘图，获得动态结果。
 
编程过程会涉及到积分方程的离散化、方程组的求解、求解结果的绘图等问题，编程和调试会花费一定的时间，建议在编写和调试Matlab程序的时候，先处理简单的结果，等程序调试通过后，再计算复杂的结果！
 
但我在离散化时还是不知怎么入手，大家帮忙想哈

function y=Broy(x0,eps)
if nargin==1
    eps=1e-6;
else if nargin<1
        error
        return
    end
x0=[50,50,50];
A=eye(length(x0));
x1=x0-f(x0)/A;
n=1;
while(norm(x1-x0)>=1e-6&(n<=10000))
x0=x1;
x1=x0-f(x0)/A;
p=x1-x0;
q=f(x1)-f(x0);
A=A+(q-p*A)'*p/norm(p);
n=n+1;
end
y=x1;
end


function t=f(x)
syms R W ;
l=(0.021*W-0.021*x(1)-0.0008*x(2)-0.021*x(3)+0.02*x(3)*R+10.9)/R*exp(-12.5*(log(R)-9/125)^2);
g=(R-1.0513)*(0.021*W-0.021*x(1)-0.0008*x(2)-0.021*x(3)+0.02*x(3)*R+10.9)/R*exp(-12.5*(log(R)-9/125)^2);
p=(-0.00084*W+0.00084*x(1)+0.00181*x(2)+0.00084*x(3)-0.0008*x(3)*R-0.4)/R*exp(-12.5*(log(R)-9/125)^2);
r=linspace(0.01,20,51);
h=(20-0.01)/50,
a(1)=subs(l,R,r(1));b(1)=subs(g,R,r(1));c(1)=subs(p,R,r(1));
a(51)=subs(l,R,r(51));b(51)=subs(g,R,r(51));c(51)=subs(p,R,r(51));
t(1)=(a(1)+a(51))/2*h;t(2)=(b(1)+b(51))/2*h;t(3)=(c(1)+c(51))/2*h;
for i=2:50
    a(i)=subs(l,R,r(i));b(i)=subs(g,R,r(i));c(i)=subs(p,R,r(i));
    t(1)=t(1)+a(i)*h;t(2)=t(2)+a(i)*h;t(3)=t(3)+a(i)*h;
end
t(1)=0.02*x(1)+10-t(1);
t=[t(1) t(2) t(3)];
end

使用这两个M文件，但主程序怎么编都有问题，希望大家指导下！