f.m文件：
% du/dx=ffunction d=f(x, u)d=(2*u/x) +(x^2)*exp(x); end
欧拉法
Euler.m文件：
function [x,u]=Euler(h)x1=1;x2=2;u1=0;N=abs(x2-x1)/h; u=zeros(N+1,1);x=zeros(N+1,1);u(1)=u1;x(1)=x1; for i=1: Nx(i+1)=x1+i*h;u(i+1)=u(i)+h*f( x(i), u(i) ); end end
控制台输入：
[x1,y1] =Euler(0.1);[x2,y2] =Euler(0.05);[x3,y3] =Euler(0.01);plot(x1,y1,'b');hold on;plot(x2,y2,'g');hold on;plot(x3,y3,'r');xlabel('x');ylabel('u');legend('h=0.1','h=0.05','h=0.01','Location','northwest');
改进的Euler法
improvedEuler.m文件：
function [x,u]=improvedEuler(h)x1=1;x2=2;u1=0;N=abs(x2-x1)/h; u=zeros(N+1,1);x=zeros(N+1,1);u(1)=u1;x(1)=x1; for i=1: Nx(i+1)=x1+i*h;up=u(i)+h*f( x(i), u(i) );uc=u(i)+h*f( x(i+1), up);u(i+1)=(up+uc)/2; end end
控制台输入：
[x1,y1] = improvedEuler(0.1);[x2,y2] = improvedEuler(0.05);[x3,y3] = improvedEuler(0.01);plot(x1,y1,'b');hold on;plot(x2,y2,'g');hold on;plot(x3,y3,'r');xlabel('x');ylabel('u');legend('h=0.1','h=0.05','h=0.01','Location','northwest');
四阶龙格库塔法
RungeKutta.m文件：
function [x,u]=RungeKutta(h)a=1;b=2;u1=0;N=floor((b-a)/h);x(1)=a;u(1)=u1;for ii=1:Nx(ii+1)=x(ii)+h;k1=f(x(ii),u(ii));k2=f(x(ii)+h/2,u(ii)+h*k1/2);k3=f(x(ii)+h/2,u(ii)+h*k2/2);k4=f(x(ii)+h,u(ii)+h*k3);u(ii+1)=u(ii)+h*(k1+2*k2+2*k3+k4)/6;end
控制台输入：
[x1,y1] = RungeKutta(0.1);[x2,y2] = RungeKutta(0.05);[x3,y3] = RungeKutta(0.01);plot(x1,y1,'b');hold on;plot(x2,y2,'g');hold on;plot(x3,y3,'r');xlabel('x');ylabel('u');legend('h=0.1','h=0.05','h=0.01','Location','northwest');
