线性多步法

clear all;
clc;
x0 = 0;
b = 1.5; 
y0 = 3;

h = 0.1;
df = @diff;
[k,X,Y,wucha,P] = Adams4x(df,x0,b,y0,h);
temp = exact(X);
plot(X,Y,'r',X,temp,'*');
%plot(X,temp,'*');
function result = exact(x)
     result = 3./(1+x.^3);%exp(-5*x);
end
function result = diff(x,y)
    result = -x.^2.*y.^2;%-5*x;%x.^3 - y./x;
end
function [k,X,Y,wucha,P] = Adams4x(funcfcn,x0,b,y0,h)
x = x0;y = y0;p = 128; n = fix((b-x0)/h);
if n<5
	return;
end
X = zeros(p,1);
Y = zeros(p,length(y)); f = zeros(p,1);
k = 1; X(k) = x; Y(k,:) = y';
for k =2:4
	c1 = 1/6; c2 = 2/6;  c3 = 2/6;
	c4 = 1/6; a2 = 1/2;  a3 = 1/2;
	a4 = 1;   b21 = 1/2; b31 = 0;
	b32 = 1/2;b41 = 0;   b42 = 0;  b43 = 1;
	x1 = x+a2*h; x2 = x+a3*h;
	x3 = x+a4*h; 
    k1 = feval(funcfcn,x,y);
	y1 = y+b21*h*k1; 
    x = x+h;
	k2 = feval(funcfcn,x1,y1);
	y2 = y+b31*h*k1+b32*h*k2;
	k3 = feval(funcfcn,x2,y2);
    y3 = y+b41*h*k1+b42*h*k2+b43*h*k3;%y3 = y+h*k3;    
	k4 = feval(funcfcn,x3,y3);
	y = y+h*(c1*k1+c2*k2+c3*k3+c4*k4); %y = y+h*(1/6*k1+...+1/6*k4);
    X(k) = x; Y(k,:) = y;
end
X;Y;f(1:4) = feval(funcfcn,X(1:4),Y(1:4));
count_max = 100;%隐函数的最大迭代次数

for k = 4:n
    X(k+1) = X(1) + h*k;
    f(k+1) = feval(funcfcn,X(k+1),Y(k));
    count = 1;
    du = 1;
    while du>1e-6 && count <count_max
        s1 = Y(k) +(h/24)*sum(((f(k-2:k+1))'.*[1 -5 19 9]));
        f(k+1) = feval(funcfcn,X(k+1),s1);
        s2 = Y(k) +(h/24)*sum(((f(k-2:k+1))'.*[1 -5 19 9]));
        Y(k+1) = s2;
        du = s2-s1;
        count = count+1;
    end
end
  %Y(k+1) = Y(k) +(h/24)*sum(((f(k-2:k+1))'.*[1 -5 19 9]));
% for k = 4:n
% 	f(k) = feval(funcfcn,X(k),Y(k));
% 	X(k+1) = X(1) + h*k;
% 	Y(k+1) = Y(k) +(h/24)*sum(((f(k-3:k))'.*[-9 37 -59 55]));
%     
% % 	f(k+1) = feval(funcfcn,X(k+1),Y(k+1));%预估校正公式
% %     Y(k+1) = Y(k) +(h/24)*sum(((f(k-2:k+1))'.*[1 -5 19 9]));
%  
% end
for k = 2:n+1
	wucha(k) = norm(Y(k) - Y(k-1));% k = k+1;
end
X = X(1:n+1);
Y = Y(1:n+1);
n = 1:n+1;
% wucha = wucha(1:n,:);
P = [n',X,Y,wucha'];
end