clear all
close all
clc
%比较Euler法与RK算法
format long

m=0.4;%input('请输入要求自变量点');
h=0.1;%input('请输入步长h=');
f = @(x,y)(8-3*y);  % 定义 f
x=[0:h:m];
n=round(m/h);
y=zeros(1,n+1);
y(1)=2; %初始值

%4nd RK
tic
for i=1:n
    K1 =f(x(i),  y(i)); 
    K2 =f(x(i) + h/2, y(i) + h*K1/2 ); 
    K3 =f(x(i) + h/2, y(i) + h*K2/2 ); 
    K4 =f(x(i) + h, y(i) + h*K3 ); 
    y(i+1)=y(i) + 1/6*h*(K1 + 2*K2 + 2*K3 + K4 )
end
toc
%Euler法
%xx=[0:h:m];
yy=zeros(1,n+1);
yy(1)=2; %初始值ocic
tic
for i=1:n
  yy(i+1)=yy(i)+h*(8-3*yy(i));
end
toc
%%精确解
%xxx=0:h:m;
yyy=@(x)8/3-2/3*exp(-3*x);

yyy(m)
error1=y(n+1)-yyy(m) %error of RK
error2=yy(n+1)-yyy(m) % Euler
figure;plot(x,y,'-r',x,yy,'-k',x,yyy(x),'b')
legend('RK方法','Euler法','精确解')
title('经典四阶RK法与Euler法比较')
xlabel('x')
ylabel('y')
figure;plot(x,y-yyy(x),'-r')%,x,yy-yyy(xxx),'-b')
title('经典四阶RK法的误差')
xlabel('x')
ylabel('error')



