matlab：使用4阶龙格库塔方法求微分方程组的值

参考书目：《数值分析（matlab版）》  作者：周璐 翻译

%调用龙格库塔方法求解微分方程组，P399,9.7 微分方程组
function [t,z] = my_rk4b(f, t0, tf , Za, h)
%f-函数； t0，tf:区间； Za，初值,是行向量；h 步长
%输出的z为一个矩阵，行数为离散点的个数+1，列表示向量

M = floor((tf - t0)/h);     % 离散点的个数
t = zeros(M +1, 1);
z = zeros(M +1, size(Za,2));        %size(Za,1)表示变量个数，每一列是个变量

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

function s = myfun_b(t,z)
%z是个向量,1*2
%输出s也是向量,1*2
     s1 = z(1,1) + 2*z(1,2);
     s2 = 3*z(1,1) + 2*z(1,2);
  s = [s1 s2];
 end

clear all;clc
%调用龙格库塔方法求