龙格库塔方法matlab实现~functionff=rk(yy,x0,y0,h,a,b)%yy为y的导函数,x0,y0,为初值,h为步长,a,b为区间c=(b-a)/h+1;i1=1;%c为迭代步数;i1为迭代步数累加值y=y0;z=zeros(c,6);%z生成c行,5列的零矩阵存放结果;%每行存放c次迭代结果,每列分别存放k1~k4及y的结果forx=a:h:bifi1=ck1=feval(yy,x,y);k2=feval(yy,x+h/2,y+(h*k1)/2);k3=feval(yy,x+h/2,y+(h*k2)/2);k4=feval(yy,x+h,y+h*k3);y=y+(h/6)*(k1+2*k2+2*k3+k4);z(i1,1)=x;z(i1,2)=k1;z(i1,3)=k2;z(i1,4)=k3;z(i1,5)=k4;z(i1,6)=y;i1=i1+1;endendfprintf(‘结果矩阵,第一列为x(n),第二列~第五列为k1~k4,第六列为y(n+1)的结果')z%在命令框输入下列语句%yy=inline('x+y');%rk(yy,0,1,0.2,0,1)%将得到结果%结果矩阵,第一列为x(n),第二列~第五列为k1~k4第六列为y(n+1)的结果%z=%01.00001.20001.22001.44401.2428%0.20001.44281.68711.71151.98511.5836%0.40001.98362.28202.31182.64602.0442%0.60002.64423.00863.04513.45322.6510%0.80003.45103.89613.94074.43923.4365%1.00004.43654.98025.03455.64344.4401