runge_kutta.m

常微分方程数值解法Runge_Kutta解法

Runge_Kutta.m
function [x,y]=Runge_Kutta(fun,x0,y0,h,N)
x=zeros(1,N+1);y=zeros(1,N+1);
x(1)=x0;y(1)=y0;
for n=1:N
    x(n+1)=x(n)+h;
    k1=feval(fun,x(n),y(n));
    k2=feval(fun,x(n)+h/2,y(n)+k1*h/2);
    k3=feval(fun,x(n)+h/2,y(n)+k2*h/2);
    k4=feval(fun,x(n)+h,y(n)+k3*h);
    y(n+1)=y(n)+1/6*h*(k1+2*(k2+k3)+k4);
end