【求助】用龙格库塔法数值计算（毕业论文相关，急。导师不管事只能网上求助了）

问题如下：

x_1=x,x_2=x’,

相应初始条件为：x_1 (0)=1,x_2 (0)=0

求：
突加载荷时，也就是P（t）=p时
（1） t-x曲线；
（2） 不同P/T/U下的x-x’曲线；
（3） 不同U/T下P-Xmax曲线

程序如下：

主程序：
function [x,y]=runge_kutta1(ufunc,y0,h2,a,b) %参数表顺序依次是微分方程组的函数名称，初始值向量，步长，时间起点，时间终点（参数形式参考了ode45函数）
n=floor((b-a)/h2); %求步数
x(1)=a; %时间起点
y(:,1)=y0; %赋初值，可以是向量，但是要注意维数
for ii=1:n
x(ii+1)=x(ii)+h2;
k1=ufunc(x(ii),y(:,ii));
k2=ufunc(x(ii)+h2/2,y(:,ii)+h2k1/2);
k3=ufunc(x(ii)+h2/2,y(:,ii)+h2k2/2);
k4=ufunc(x(ii)+h2,y(:,ii)+h2k3);
y(:,ii+1)=y(:,ii)+h2(k1+2k2+2k3+k4)/6; %按照龙格库塔方法进行数值求解
end
end

function [C,L]=lagran(X,Y)
% input - X is a vector that contains a list of abscissas
% - Y is a vector that contains a list of ordinates
% output - C is a matrix that contains the coefficients of the lagrange interpolatory polynomial
%- L is a matrix that contains the lagrange coefficients polynomial
w=length(X);
n=w-1;
L=zeros(w,w);
for k=1:n+1
V=1;
for j=1:n+1
if k~=j
V=conv(V,poly(X(j)))/(X(k)-X(j));%对多项式做卷积运算
end
end
L(k,:)=V;
end
C=Y*L;
end

function [X,Y]=mid(q,A,B,h)
i=0;
for x=1:h:q
F_sum=0;i=i+1;X(i)=x;
N=100;
miu=[1.74e6,3.501e3,-2.684e4];a=[1.3,5,-2];
delta=(B/A)^3-1;
limmin=x;limmax=((delta+x^{3)/(delta+1))}(1/3);%积分上下限
xi=linspace(limmin,limmax,N); %linspace(x1,x2,N)用于产生x1,x2之间的N点行线性的矢源量，其中x1、x2、N分别为起始值、终止值、元素个数。百若默认N，默认点数为100.
dxi=(limmax-limmin)/(N-1);
expr=zeros(1,N);
for p=1:3
expr=expr+miu§.*(xi.^(-a§)-xi.(1/2.a§)); %.(xi.^{(-a§)*lambda_z}(-a§)-xi.^a§);
end
expr=expr./xi./2;
f_sum=sum(expr)*dxi;
F_sum=F_sum+f_sum;
Y(i)=F_sum;
end
end

function dx=test_fun1(t,x,flag,P,U,wen,A,B)
R=0.039;
%P=50000;
c4=-269.982;
% wen=300;
wen0=300;
delta=(B/A)^3-1;
rho=950;epsilong=2.21e-11;
[X,Y]=mid(1.10,0.038,0.04,0.02);
S=polyfit(X,Y,2);%使用polyfit()函数获得通过当前数据集的拟合结果，最后一个参数为多项式的最高次数。

%之后用polyval()函数获得当前数据集对应的输出结果。
f=0;
for i=1:3
f=S(i)*x(1)^(3-i)+f;
end

dail= delta+x(1)^3;
s=U^{2*epsilong*x(1)}2*(dail)^(2/3)(( dail)^(-1)-x(1)(-3))/3/A/((( dail)^(1/3)-x(1)))2;
k=4/3c4*(wen-wen0)(R^4*A(-3)(dail^(-1)-x(1)(-3))-R*(-2)A^3delta);
%P=p(t);
F=(-f-P+s+k)/rho/A^2;
J=x(2)^2/2*(x(1)4*dail^{(-4/3)-4*x(1)*dail}(1/3)+3);
dx = zeros(2,1); %初始化列向量
dx(1)=x(2);
dx(2)=(J+F)/x(1)^2;

end

t-x曲线程序：

clc;clear;
P=100000;
U=10000;
wen=320;A=0.038;B=0.040;
[T,Y]=ode45(‘test_fun1’,[0 0.1],[1,0],[],P,U,wen,A,B);
plot(T,Y(:,1),’-’)
xlabel(‘Time:t/s’,‘FontSize’,20,‘FontWeight’,‘bold’);
ylabel(‘Deformation:x’,‘FontSize’,20,‘FontWeight’,‘bold’);
title(‘时程曲线’,‘FontSize’,20,‘FontWeight’,‘bold’,‘FontName’,‘楷书’)
hleg1=legend([‘P=’,num2str§,‘Pa,U=’,num2str(U),‘V,T=’,num2str(wen),‘K’]);

想问问以上程序哪里出了问题呢？还有p-Xmax曲线的程序实在不会写，怎么算都是错的，求大神教一教，是毕业论文要用的，比较急，谢谢各位了。