rk.m

来自「matlab codes for solving ODEs, including」· M 代码 · 共 18 行

function [tout,yout]=rk(f,tspan,y0,N)

h=(tspan(2)-tspan(1))/N;
t = tspan(1);
tout = t;
y=y0(:);
yout = y.';
for n=1:N
    k1 = feval(f,t, y);
    k2 = feval(f,t+0.5*h, y+0.5*h*k1);
    k3 = feval(f,t+0.5*h, y+0.5*h*k2);
    k4 = feval(f,t+h, y+h*k3);
    y = y+h*(k1+2*k2+2*k3+k4)/6;
    t = t+h;
    yout=[yout;y.'];
    tout=[tout;t];
end