第八章常微方程的初值问题1.pptVIP

function [x,y]=RK4(fun,x0,y0,xmax,h) %fun为显式一阶微分方程右端的函数 %x0,y0为初始条件,满足y(x0)=y0 %xmax为x的取值最大值 %h为将x等分后的步长 N=(xmax-x0)/h+1;%N为总的节点数 x(1)=x0;y(1)=y0; for n=1:N-1 x(n+1)=x(n)+h; k1=h*feval(fun,x(n),y(n)); k2=h*feval(fun,x(n)+1/2*h,y(n)+k1/2); k3=h*feval(fun,x(n)+1/2*h,y(n)+k2/2); k4=h*feval(fun,x(n+1),y(n)+k3); y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4); end fun=inline(y+2*x/y^2,x,y); [x1,y1]=Euler_mend(fun,0,1,2,0.1); [x2,y2]=RK4(fun,0,1,2,0.1); x=0:0.1:2; y=1/3*(-18-54*x+45*exp(3*x)).^(1/3); plot(x1,y1,*,x2,y2,s,x,y) axis([0,2,0,9]) 再分别调用Eul