范德堡方程的龙格库塔算法求解

龙格库塔算法：function[x,y]=lgkt(f,x0,xf,y0,h)h=0.1;n=fix((xf-x0)/h);x(1)=x0;y(:,1)=y0;fori=1:n-1x(i+1)=x(i)+h;k1=f(x(i),y(:,i));k2=f(x(i)+h/2,y(:,i)+h*k1/2);k3=f(x(i)+h/2,y(:,i)+h*k2/2);k4=f(x(i)+h,y(:,i)+h*k3);y(:,i+1)=y(:,i)+h*(k1+2*k2+2*k3+k4)/6;end函数化简：functionf=f(t,y)f=[y(2);(1-y(1).^2).*y(2)-y(1)];CommandWindow的输入：x0=input('x0=');xf=input('xf=');whilex0=xfdisp('输入不合逻辑，请重新输入')x0=input('x0=');xf=input('xf=');endy0=[18];h=0.1;[t,F1]=lgkt(@f,x0,xf,y0,h)plot(t,F1)gridtitle('范德堡的龙格库塔算法之求解')legend('y原函数，y的一阶导')