解常微分方程的多步法（有admas外插法，admas内插法，一般多步法（Milne法））

解常微分方程的多步法（有admas外插法，admas内插法，一般多步法（Milne法））

admas外插法

%admas外插法，具有k-1阶的精度,下面编写k等于3的外插公式
clear
clc
syms x y
f=3*x^2;
a=0;
b=10;
y0=0;
y1=0.001;
y2=0.008;
y3=0.027;
h=0.1;
n=(b-a)/h+1;
xx=zeros(n,1);
yy=zeros(n,1);
yy=zeros(n,1);
yy(1)=y0;
yy(2)=y1;
yy(3)=y2;
yy(4)=y3;
for i=1:n
    xx(i)=a+h*(i-1); 
end
for i=5:n
   k1=55/24*h*(subs(f,{x,y},{xx(i-1),yy(i-1)}));
   k2=-59/24*h*(subs(f,{x,y},{xx(i-2),yy(i-2)})); 
   k3=37/24*h*(subs(f,{x,y},{xx(i-3),yy(i-3)}));
   k4=-9/24*h*(subs(f,{x,y},{xx(i-4),yy(i-4)}));
   yy(i)=k1+k2+k3+k4+yy(i-1);
end
ff1=[xx,yy];

admas内插法

%admas内插法，具有k阶的精度,但是为隐式方程，所以需要求解方程，下面编写当p等于1，k=4的admas内插法
clear
clc
syms x y
f=3*x^2;
a=0;
b=10;
y0=0;
y1=0.001;
y2=0.008;
y3=0.027;
h=0.1;
n=(b-a)/h+1;
xx=zeros(n,1);
yy=zeros(n,1);
yy=zeros(n,1);
yy(1)=y0;
yy(2)=y1;
yy(3)=y2;
yy(4)=y3;
for i=1:n
    xx(i)=a+h*(i-1); 
end
for i=5:n
   syms p
   k1=subs(f,{x,y},{xx(i),p})*h*251;
   k2=subs(f,{x,y},{xx(i-1),yy(i-1)})*h*646;
   k3=subs(f,{x,y},{xx(i-2),yy(i-2)})*h*(-264);
   k4=subs(f,{x,y},{xx(i-3),yy(i-3)})*h*106;
   k5=subs(f,{x,y},{xx(i-4),yy(i-4)})*h*(-19);
   yy(i)=double(solve(p==yy(i-1)+k1+k2+k3+k4+k5,p));
end
ff1=[xx,yy];

一般多步法（Milne法）

%一般线性多步法，下面介绍Milne法
clear
clc
syms x y
f=3*x^2;
a=0;
b=10;
y0=0;
y1=0.001;
y2=0.008;
y3=0.027;
h=0.1;
n=(b-a)/h+1;
xx=zeros(n,1);
yy=zeros(n,1);
yy=zeros(n,1);
yy(1)=y0;
yy(2)=y1;
for i=1:n
    xx(i)=a+h*(i-1); 
end
for i=3:n
   syms p
   k1=h/3*(subs(f,{x,y},{xx(i),p})+4*subs(f,{x,y},{xx(i-1),yy(i-1)})+subs(f,{x,y},{xx(i-2),yy(i-2)}));
   yy(i)=double(solve(p==yy(i-2)+k1,p));
end
ff1=[xx,yy];