gusucode.com > 高等数学问题求解源码程序 > CH14/Classical_RK4.m

    function [x,y]=Classical_RK4(odefun,xspan,y0,h,varargin)
%CLASSICAL_RK4   经典四阶Runge-Kutta法求解初值问题的数值解
% [X,Y]=CLASSICAL_RK4(ODEFUN,XSPAN,Y0,H)  经典四阶Runge-Kutta法求解微分方程ODEFUN
% [X,Y]=CLASSICAL_RK4(ODEFUN,XSPAN,Y0,H,P1,P2,...)  经典四阶Runge-Kutta法求解
%                                  微分方程ODEFUN，ODEFUN包含附加参数P1,P2,...
%
% 输入参数：
%     ---ODEFUN：微分方程的函数描述
%     ---XSPAN：求解区间[x0,xn]
%     ---Y0：初始条件
%     ---H：迭代步长
%     ---P1,P2,...：ODEFUN函数的附加参数
% 输出参数：
%     ---X：返回的节点，即X=XSPAN(1):H:XSPAN(2)
%     ---Y：微分方程的数值解
%
% See also Explicit_Euler

x=xspan(1):h:xspan(2);
N=length(x);
y=zeros(1,N);
y(1)=y0;
for k=1:length(x)-1
    K1=feval(odefun,x(k),y(k),varargin{:});
    K2=feval(odefun,x(k)+h/2,y(k)+h/2*K1,varargin{:});
    K3=feval(odefun,x(k)+h/2,y(k)+h/2*K2,varargin{:});
    K4=feval(odefun,x(k)+h,y(k)+h*K3,varargin{:});
    y(k+1)=y(k)+h/6*(K1+2*K2+2*K3+K4);
end
x=x';y=y';
web -broswer http://www.ilovematlab.cn/forum-221-1.html