% 例题3-1：求微分方程:Euler法、泰勒级数法、改进的Euler法
clear
clc;close all
format long

x_in=2; %要求点的自变量值
y=@(x)exp(-0.5*x.^2); %精确值
x=0:0.001:x_in;

h=0.1;
x0=0;
m=(x_in-x0)/h;
xx=x0:h:m*h+x0;
%--------------改进的Euler法----------------
y_2rk=zeros(1,m+1);
y_2rk(1)=1;
f=@(x,y)(-x*y);
%while 
for n=1:m
    k1=f(xx(n),y_2rk(n));
    k2=f(xx(n)+h,y_2rk(n)+h*k1);
    y_2rk(n+1)=y_2rk(n)+h/2*(k1+k2);
    error=abs(y_2rk(n+1)-(y_2rk(n)+h*k1));
    while error>10^-7
        k1=f(xx(n),y_2rk(n));
        k2=f(xx(n)+h,y_2rk(n+1));
        temp=y_2rk(n+1);
        y_2rk(n+1)=y_2rk(n)+h/2*(k1+k2);
        error=abs(y_2rk(n+1)-temp);
    end
end

figure;plot(x,y(x),'r',xx,y_2rk,'b') 
legend('精确解','自洽循环的改进的Euler法')
xlabel('x')
ylabel('y')
title('采用自洽循环的改进的Euler法求解常微分方程的初值问题')


figure; plot(xx,y(xx)-y_2rk,'r')
legend('自洽循环的改进的Euler法')
xlabel('x')
ylabel('error')
title('自洽循环的改进的Euler法的误差')

fprintf('y(x_in)的精确值为：%f\n',y(x_in));
fprintf('改进的Euler法求得的值为：%f\n',y_2rk(m+1));
error2=y(x_in)-y_2rk(m+1);
fprintf('改进的Euler法求得的值为：%f\n',error2);

