原方程、文献图如图所示。其中R=1;V=1;a=0.1;我模拟的微分模拟图不知道为啥是这样的。
微分方程组：

function du=f(xi,u)
%a替换alpha
syms R V xi a u0 u1 u2 u3
du=zeros(4,1);
du0=-u0-1./2*V*u1;   %22 初始值
du1=-V*u0-(1+R*exp(-a*xi)*cosh(a*V *xi))*u1-(1./2*V-R*exp(-a*xi)*sinh(a*V *xi))*u2;   %23;a替换alpha
du2=-(1./2*V-R*exp(-a*xi)*sinh(a*V *xi))*u1-(1+4*R*exp(-a*xi)*cosh(a*V *xi))*u2-(1./2*V-3*R*exp(-a*xi)*sinh(a*V *xi))*u3;   %24
du3=-(1./2*V-R*exp(-a*xi)*sinh(a*V *xi))*u2-(1+9*R*exp(-a*xi)*cosh(a*V *xi))*u3';   %25
end

R=1;
V=1;
a=0.1;
[xi,u]=ode45('f',[0 10],[100 0 0 0]);   %求数值解
plot(xi,u(:,1),'-r',xi,u(:,2),'.k',xi,u(:,3),'*g',xi,u(:,4),':b')

积分方程组是基于上述结果的，也不知怎么搞：
function N=integral(xi)
syms xi V N u0 u1 u2 u3 N0 N1 N2 N3
N0=u0+1./2*V*u1;
N1=V*u0+u1+1./2*V*u2;
N2=1./2*V*u1+u2+1./2*V*u3;
N3=1./2*V*u2+u3
end

N0=quad(@integral,0,10);
N1=quad(@integral,0,10);
N2=quad(@integral,0,10);
N3=quad(@integral,0,10);
plot(xi,N(:,1),'-r',xi,N(:,2),'.k',xi,N(:,3),'*g',xi,N(:,4),':b')
不知怎么发不了图片，图片：https://zhidao.baidu.com/questio ... ?entry=qb_uhome_tag 返回小木虫查看更多