t(i+1)=t(i)+h;
y(i+1)=y(i)+1/2*h*(fun(t(i),y(i))+ fun(t(i+1),y1(i+1)))
end
plot(t,y,'*r')
function y=fun(t,y);
y=y+1;
调用：Euler(0,3,[0,2],
得到解析解：hold on;
y=dsolve('Dy=y+1','(y(0)=3)','t');
t=t0;
x(1)=int(1);
y(1)=int(2);
for i=1:n
x1(i+1)=x(i)+h*xfun(t(i),x(i),y(i));
y1(i+1)=y(i)+h*yfun(t(i),x(i),y(i));
t(i+1)=t(i)+h;
x(i+1)=x(i)+1/2*h*(xfun(t(i),x(i),y(i))+xfun(t(i+1),x1(i+1),y1(i+1)));
y(i+1)=y(i)+1/2*h*(yfun(t(i),x(i),y(i))+yfun(t(i+1),x1(i+1),y1(i+1)));
end
plot(t,x,'o-r')
hold on
plot(t,y,'*-g&#
function x=xfun(t,x,y);
x=y+y^2;
ezplot(y,[0,2])
图像：
代码：
function Euler1(t0,y0,inv,h)
n=round(inv(2)-inv(1))/h;
t(1)=t0;
y(1)=y0;
for i=1:n
y1(i+1)=y(i)+h*fun(t(i),y(i));
t(i+1)=t(i)+h;
y(i+1)=y(i)+1/2*h*(fun(t(i),y(i))+ fun(t(i+1),y1(i+1)))
hold on
plot(t,y,'*-g')
hold on
plot(x,y)
function x=xfun(t,x,y);
x=y;
function y=yfun(t,x,y);
y=-2*x-3*y;
调用函数：Euler4(0,[1,1],5,
图像：
（2）代码：function Euler5(t0,int,n,h)
end
plot(t,y,'*r')
function y=fun(t,y);
y=y^2-4*t;
调用：
Euler1(0,,[0,2],
图像：
代码：
function Euler2(t0,y0,inv,h)
n=round(inv(2)-inv(1))/h;
t(1)=t0;
y(1)=y0;
for i=1:n
y1(i+1)=y(i)+h*fun(t(i),y(i));
t(i+1)=t(i)+h;
y(i+1)=y(i)+1/2*h*(fun(t(i),y(i))+ fun(t(i+1),y1(i+1)))
end
plot(t,y,'*r')
function y=fun(t,y);
y=(3-y)*(y+1);
function y=yfun(t,x,y);
y=-x+*y-x*y+*y^2;
调用函数：Euler5(0,[1,1],5,
图像：
调用：
Euler2(0,4,[0,5],1)
得到解析解：
hold on;
y=dsolve('Dy=(3-y)*(y+1)','y(0)=4','t');
ezplot(y)
图像：
代码：
function Euler2(t0,y0,inv,h)
n=round(inv(2)-inv(1))/h;
t(1)=t0;
（1）代码：
function Euler4(t0,int,n,h)
t=t0;
x(1)=int(1);
y(1)=int(2);
for i=1:n
x1(i+1)=x(i)+h*xfun(t(i),x(i),y(i));
y1(i+1)=y(i)+h*yfun(t(i),x(i),y(i));
t(i+1)=t(i)+h;
y=(3-y)*(y+1);
调用：
Euler2(0,4,[0,5],
得到解析解：
hold on;
y=dsolve('Dy=(3-y)*(y+1)','y(0)=4','t');
ezplot(y)
图像：
14.考虑满足初始条件（x(0)，y(0)）=(1,1)的下列方程组：
选定时间步长 t=，n=5.用改进欧拉方法求两个方程组的近似解；
P77 31.利用改进欧拉方法计算下列初值问题，并画出近似解的草图：
代码：
%改进欧拉法
function Euler(t0,y0,inv,h)
n=round(inv(2)-inv(1))/h;
t(1)=t0;
y(1)=y0;
for i=1:n
y1(i+1)=y(i)+h*fun(t(i),y(i));
y(1)=y0;
for i=1:n
y1(i+1)=y(i)+h*fun(t(i),y(i));
t(i+1)=t(i)+h;
y(i+1)=y(i)+1/2*h*(fun(t(i),y(i))+ fun(t(i+1),y1(i+1)))
end
plot(t,y,'*r')
function y=fun(t,y);
x(i+1)=x(i)+1/2*h*(xfun(t(i),x(i),y(i))+xfun(t(i+1),x1(i+1),y1(i+1)));
y(i+1)=y(i)+1/2*h*(yfun(t(i),x(i),y(i))+yfun(t(i+1),x1(i+1),y1(i+1)));
end
plot(t,x,'o-r')