经典R-K

yp.m:

function dxdy=yp(x,y)
dxdy=-2*y-4*x;

ExpEulMtd.m:

%Explicit Euler method
%ExpEulMtd

h=0.1;
x=0:h:0.5;
[row,col]=size(x);
y=zeros(row,col);
i=1;
y(i)=2;
for i=2:col
    y(i)=y(i-1)+h*yp(x(i-1),y(i-1));
end
figure
plot(x,y)
y

ClassicRKMtd.m

%Classic R-K method
%ClassicRKMtd
h=0.1;
x=0:h:0.5;
[row,col]=size(x);
y=zeros(row,col);
i=1;
y(i)=2;
for i=2:col
    k1=yp(x(i-1),y(i-1));    
    y_mid=y(i-1)+h*k1/2;
    x_mid=x(i-1)+h/2;
    k2=yp(x_mid,y_mid); 
    y_mid=y(i-1)+h*k2/2;
    k3=yp(x_mid,y_mid); 
    y(i)=y(i-1)+h*k3;
    k4=yp(x(i),y(i));
    k_avg=(k1+2*k2+2*k3+k4)/6;
    y(i)=y(i-1)+h*k_avg;    
end
figure
plot(x,y)
y

trueSovle.m

%trueSovle
h=0.1;
x=0:h:0.5;
[row,col]=size(x);
y=zeros(row,col);
y=exp(-2*x)-2*x+1;
figure
plot(x,y)
y

>> ClassicRKMtd


y =


    2.0000    1.6187    1.2703    0.9488    0.6493    0.3679


>> trueSovle


y =


    2.0000    1.6187    1.2703    0.9488    0.6493    0.3679


>>