四阶龙格库塔法解数值微分.zip

clc clear s=input('请按格式输入f(x,y)=','s'); f=inline(s,'x','y'); a=input('请按格式输入x的取值下限a='); b=input('请按格式输入x的取值上限b='); ya=input('请按格式输入ya的值ya='); h=input('请按格式输入步长h='); n=(b-a)/h; y=ones(1,n+1); x=ones(1,n+1); y(1)=ya; for i=1:1:n+1 x(i)=a+(i-1)*h; end for i=1:1:n k1=h*f(x(i),y(i)); k2=h*f(x(i)+0.5*h,y(i)+0.5*k1); k3=h*f(x(i)+0.5*h,y(i)+0.5*k2); k4=h*f(x(i)+h,y(i)+k3); y(i+1)=y(i)+1/6*(k1+2*k2+2*k3+k4); fprintf('y(%d)=%.4f \n',i,y(i+1)); end