хотелрешить систему из трех дифференциальных уравнений методом Рунге Кутты 4 в Matlab (Ode45 не допускается).

После долгих поисков я смог найти в Интернете либо неразборчивые примеры, либо общие объяснения, которые вообще не содержат примеров. Я хотел бы получить конкретный пример того, как правильно реализовать свое решение, или решение сопоставимой проблемы, на которой я могу основываться.

Я пришел довольно далеко;мой текущий код выплевывает матрицу с 2 правильными десятичными знаками на большинстве компонентовкоторым я вполне доволен.

Однако, когда размер шага уменьшается, ошибки становятся огромными. Я знаю, что созданный цикл не совсем корректен. Возможно, я неправильно определил функции, но я вполне уверен, что проблема будет решена, если будут внесены некоторые незначительные изменения в цикл for, потому что он, похоже, довольно хорошо решает систему уравнений уже в своем текущем состоянии.

clear all, close all, clc

%{
____________________TASK:______________________
Solve the system of differential equations below 
in the interval 0<t<1, with stepsize h = 0.1.
x'= y                 x(0)=1
y'= -x-2e^t+1         y(0)=0   ,   where x=x(t), y=y(t),  z=z(t)  
z'= -x - e^t + 1      z(0)=1

THE EXACT SOLUTIONS for x y and z can be found in this pdf:
archives.math.utk.edu/ICTCM/VOL16/C029/paper.pdf
_______________________________________________

%}

h = 0.1;
t    = 0:h:1
N = length(t);

%Defining the functions
x    = zeros(N,1);%I am not entierly sure if x y z are supposed to be defined in this way.
y    = zeros(N,1)
z    = zeros(N,1)

f = @(t, x, y, z) -x-2*exp(t)+1;%Question: Do i need a function for x here as well??
g = @(t, x, y, z) -x - exp(t) + 1;

%Starting conditions
x(1) = 1; 
y(1) = 0;
z(1) = 1;

for i = 1:(N-1)
    K1     = h * ( y(i));%____I think z(i) is supposed to be here, but i dont know in what way. 
    L1     = h * f( t(i)          , x(i)        , y(i) ,      z(i));
    M1     = h * g( t(i)          , x(i)        , y(i) ,      z(i));

    K2     = h *  (y(i) + 1/2*L1 + 1/2*M1);%____Again, z(i) should probably be here somewhere. 
    L2     = h * f(t(i) + 1/2*h,  x(i)+1/2*K1 , y(i)+1/2*L1 , z(i)+1/2*M1);
    M2     = h * g(t(i) + 1/2*h,  x(i)+1/2*K1 , y(i)+1/2*L1 , z(i)+1/2*M1);

    K3     = h *  (y(i) + 1/2*L2 + 1/2*M2);%____z(i). Should it just be added, like "+z(i)" ? 
    L3     = h * f(t(i) + 1/2*h,  x(i) + 1/2*K2 , y(i) + 1/2*L2 , z(i) + 1/2*M2);
    M3     = h * g(t(i) + 1/2*h,  x(i) + 1/2*K2 , y(i) + 1/2*L2 , z(i) + 1/2*M2);

    K4     = h *  (y(i) + L3 + M3);%_____z(i) ... ? 
    L4     = h * f( t(i)+h    ,  x(i)+K3     , y(i)+L3,     z(i)+M3);
    M4     = h * g( t(i)+h    ,  x(i)+K3     , y(i)+L3,     z(i)+M3);

    x(i+1) = x(i)+1/6*(K1+2*K2+2*K3+K4);                                                                             
    y(i+1) = y(i)+1/6*(L1+2*L2+2*L3+L4);
    z(i+1) = z(i)+1/6*(M1+2*M2+2*M3+M4);
end

Answer_Matrix = [t' x y z]