% 四阶龙格库塔法
h=0.01; %步长
T=10;
N=T/h+1;
% 给迭代变量赋初值
Y=[0,1,0];
x=zeros(N,1); %每维N个迭代值
y=zeros(N,1);
z=zeros(N,1);
k0=zeros(1,3); %存储每次迭代产生的k0值的数组，大小为d维
k1=zeros(1,3);
k2=zeros(1,3);
k3=zeros(1,3);
k4=zeros(1,3);
i=1; %初次迭代时
% 循环
for t1=0:h:T
x(i)=Y(1);
y(i)=Y(2);
z(i)=Y(3);
k1=liu(Y); %调用liu函数
k2=liu(Y+h/2.*k1);
k3=liu(Y+h/2.*k1);
k4=liu(Y+h.*k3);
K=k1+2.*k2+2.*k3+k4; %runge-kutta法
Y=Y+h.*K.*1/6;
i=i+1;
end
% 画图
figure(4)
plot(x,t);
figure(1)
plot3(x,y,z)
figure(2);
subplot(3,1,1)
plot(x,y);
subplot(3,1,2)
plot(x,z);
subplot(3,1,3)
plot(y,z);
