f=inline('-y+x+1','x','y');%微分方程的右边项

　　dx=0.05;%x方向步长

　　xleft=0;%区域的左边界

　　xright=10;%区域的右边界

　　xx=xleft:dx:xright;%一系列离散的点

　　n=length(xx);%点的个数

　　y0=1;

　　%%(1)欧拉法

　　Euler=y0;

　　fori=2:n

　　Euler(i)=Euler(i-1)+dx*f(xx(i-1),Euler(i-1));

　　end

　　%%(2)龙格库塔法

　　RK=y0;

　　fori=2:n

　　k1=f(xx(i-1),RK(i-1));

　　k2=f(xx(i-1)+dx/2,RK(i-1)+k1*dx/2);

　　k3=f(xx(i-1)+dx/2,RK(i-1)+k2*dx/2);

　　k4=f(xx(i-1)+dx,RK(i-1)+k3*dx);

　　RK(i)=RK(i-1)+dx*(k1+2*k2+2*k3+k4)/6;

　　end

　　%%Euler和Rk法结果比较

　　plot(xx,Euler,xx,RK)

　　holdon

　　%精确解用作图

　　symsx

　　rightsolve=dsolve('Dy=-y+x+1','y(0)=1','x');%求出解析解

　　rightdata=subs(rightsolve,xx);%将xx代入解析解,得到解析解对应的数值

　　plot(xx,rightdata,'r*')

　　legend('Euler','Runge-Kutta','analytic')