yugueuler.m

利用隆格库塔法进行的计算

function[x,y]=yugueuler(fun,xspan,y0,N)
%fun(x)为初值问题的方程，xspan为区间[a,b],y0为初始值，N为划分个数
h=(xspan(2)-xspan(1))/N;
x=xspan(1):h:xspan(2);
y(1)=y0;
for n=1:length(x)-1
    K1=feval(fun,x(n),y(n));
    y(n+1)=y(n)+h*K1;
    K2=feval(fun,x(n+1),y(n+1));
    y(n+1)=y(n)+h*(K1+K2)/2;
end
x=x';y=y';