P77 31.利用改进欧拉方法计算下列初值问题，并画出近似解的草图：dy
+
=t
=
t
y
y
≤
≤
,2
;5.0
0,3
)0(
)1(=
,1
∆
dt
代码：
%改进欧拉法
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));
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],0.5)
得到解析解：hold on;
y=dsolve('Dy=y+1','(y(0)=3)','t');
ezplot(y,[0,2])
图像：
dy y
=t
-
t
y
,2
t
=
;2.0
≤
0,5.0
,4
)0(
)2(2=
≤
∆
dt
代码：
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))) end
plot(t,y,'*r')
function y=fun(t,y);
y=y^2-4*t;
调用：
Euler1(0,0.5,[0,2],0.2)
图像：
dt
代码：
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);
调用：
Euler2(0,4,[0,5],1)
得到解析解：
hold on;
y=dsolve('Dy=(3-y)*(y+1)','y(0)=4','t');
ezplot(y)
图像：
dt
代码：
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);
调用：
Euler2(0,4,[0,5],0.5)
得到解析解：
hold on;
y=dsolve('Dy=(3-y)*(y+1)','y(0)=4','t');
ezplot(y)
图像：
14.考虑满足初始条件（x(0)，y(0)）=(1,1)的下列方程组：
⎪⎪⎩⎪⎪⎨⎧+-+-=+=⎪⎪⎩⎪⎪⎨⎧--==;2.12.0,)2(;32,)1(22y xy y x dt
dy y y dt dx y x dt dy y dt dx 选定时间步长∆t=0.25，n=5.用改进欧拉方法求两个方程组的近似解；
（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;
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')
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,0.25)
图像：
（2）代码：function Euler5(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;
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')
hold on
plot(x,y)
function x=xfun(t,x,y);
x=y+y^2;
function y=yfun(t,x,y);
y=-x+0.2*y-x*y+1.2*y^2;
调用函数：Euler5(0,[1,1],5,0.25)
图像：。