四阶龙格-库塔法.zip

function [xout,yout] = RK4(odefile,xspan,y0) x0 = xspan(1); xh = xspan(2); if length(xspan) >= 3 h = xspan(3); else h = (xspan(2)-xspan(1))/100; end xout = [x0:h:xh]'; yout = []; for x = xout' K1 = eval([odefile '(x,y0)']); K2 = eval([odefile '(x+h/2,y0+0.5*K1*h)']); K3 = eval([odefile '(x+h/2,y0+0.5*K2*h)']); K4 = eval([odefile '(x+h,y0+K3*h)']); y0 = y0+(K1+2*K2+2*K3+K4)*h/6; yout=[yout; y0']; end