matlab：采用4阶龙格库塔方法求解微分方程

参考书目：《数值方法（matlab）》，作者：周璐 翻译。

%龙格-库塔方法求解微分方程
function [t,y] = my_rk4(f, t0, tf , y0, h)
%f-函数； t0，tf:区间； y0，初值；h 步长

M = floor((tf - t0)/h);     % 离散点的个数
t = zeros(M +1, 1);
y = zeros(M +1, 1);

t = [t0 : h :tf]';
y(1) = y0;
for i =1:M
    k1 = feval(f, t(i), y(i));
    k2 = feval(f, t(i)+h/2, y(i)+h/2*k1);
    k3 = feval(f, t(i)+h/2, y(i)+h/2*k2);
    k4 = feval(f, t(i)+h, y(i)+h*k3);
    y(i+1) = y(i) + h/6*( k1 + 2* k2 + 2* k3 + k4 );
end

clear all;clc
%调用龙格库塔方法求解微分方程-
%函数原型 function [t,y] = my_rk4(f, t0, tf , y0, h)
format long
t0 = 0;tf = 3;
y0 = 1;
h1 = 1; h2 = 0.5; h3 = 0.25; h4 = 1/8;

 [t1,y1] = my_rk4(@myfun, t0, tf , y0, h1)
 [t2,y2] = my_rk4(@myfun, t0, tf , y0, h2)
[t3,y3] = my_rk4(@myfun, t0, tf , y0, h3)
 [t4,y4] = my_rk4(@myfun, t0, tf , y0, h4)
 
plot(t1,y1,'b>',t2,y2,'y*' ,t3,y3,'ro',t4,y4,'g--');hold on;
 xlabel('t');ylabel('y');grid on;
 legend('h=1', 'h=1/2','h=1/4','h=1/8');
 
 plot(t4, 3*exp(-t4/2) - 2 +t4);

%第二种写法
% myfun1 = @(t,y)1/2*(t - y)
% [t,y1] = my_euler(myfun1, t0, tf , y0, h1)
% [t,y2] = my_euler(myfun1, t0, tf , y0, h2)
% [t,y3] = my_euler(myfun1, t0, tf , y0, h3)
% [t,y4] = my_euler(myfun1, t0, tf , y0, h4)

结果：