四阶Runge-Kutta方法求解常微分方程（含二阶导数）

xhx2010

以下针对著名VanderPol方程：
function s=f(t,y)
s=zeros(size(y));
s(1)=y(2);
s(2)=-y(2)*(y(1)^2-1)-y(1);

function RK_vanderpol
t0=0;tN=20;%区间
y0=[0.25;0];%初始值
h=0.001;
t=t0:h:tN;
N=length(t);
j=1;
for i=1:N
    t1=t0+h;
    k1=RK_order4(t0,y0);
    k2=RK_order4(t0+h/2,y0+h*k1/2);
    k3=RK_order4(t0+h/2,y0+h*k2/2);
    k4=RK_order4(t0+h,y0+h*k3);
    y1=y0+(h/6)*(k1+2*k2+2*k3+k4);
    yy1(j)=y1(1);
    yy2(j)=y1(2);
    t0=t1;
    y0=y1;
    j=j+1;
end

请问各位高手上面的程序怎样才能顺利的在MATLAB中运行，麻烦各位给指点下！

gentlefu

你得有RK_order4这个函数才能算。。。

xhx2010

感谢gentlefu兄弟，但是我试过了
function s=RK_order4(t,y)
s=zeros(size(y));
s(1)=y(2);
s(2)=-y(2)*(y(1)^2-1)-y(1);
弄成RK_order4.m文件，但是还是不能算，请指教

zhouyanqing1981

% function RK_vanderpol
RK_order4=@s    ;
t0=0;tN=20;%区间
y0=[0.25;0];%初始值
h=0.001;
t=t0:h:tN;
N=length(t);
j=1;
for i=1:N
    t1=t0+h;
    k1=RK_order4(t0,y0);
    k2=RK_order4(t0+h/2,y0+h*k1/2);
    k3=RK_order4(t0+h/2,y0+h*k2/2);
    k4=RK_order4(t0+h,y0+h*k3);
    y1=y0+(h/6)*(k1+2*k2+2*k3+k4);
    yy1(j)=y1(1);
    yy2(j)=y1(2);
    t0=t1;
    y0=y1;
    j=j+1;
end
%好像可以了