numerical methods - Ode-system output gives null values in Matlab

Question

I am trying to compute the numerical solution (over the interval [0,2]) for this system of differential equations using the Runga Katta (4th Order) Method and a step of my choice.

% This is not just equations, not MATLAB code
dy1/dt = 1.3 * (y2-y1) + 1.04 * 10^4 * k * y2
dy2/dt = 1.88 * 10^3 * (y4-(1+k))
dy3/dt = 1752+(266.1 * y2)-(269.3 * y3)
dy4/dt = 0.1+(320 * y2)-(321 * y4)
k = 6 * 10^(-4) * exp(20.7 - 15 * 10^3/y1)
y0 = transpose([759.167,0,600,0.1])

I adopted the following code, from this answer which ran smoothly but output a table of zeroes (barring a lone entry). I tried different step-sizes but to no avail.

clc

h = 0.2; % Step size.
t = (0:h:2)'; 
N = length(t);

y1 = zeros(N,1);    
y2 = zeros(N,1);    
y3 = zeros(N,1);
y4 = zeros(N,1);

% Initial conditions
y1(1) = 759.167;     
y2(1) = 0;    
y3(1) = 600;
y4(1) = 0.1;

% Initialise derivative functions
dy1 = @(t, y1, y2, y3, y4) (1.3*(y2-y1))+1.04*10^4*6*10^(-4)*exp(20.7 - (15*10^(3)/y1))*y2;
dy2 = @(t, y1, y2, y3, y4) 1.88*10^3*(y4-(1 + 6*10^(-4)*exp(20.7 - (15*10^(3)/y1))));
dy3 = @(t, y1, y2, y3, y4) 1752+(266.1*y2)-(269.3*y3);
dy4 = @(t, y1, y2, y3, y4) 0.1+(320*y2)-(321*y4);

% Initialise K vectors
ky1 = zeros(1,4);
ky2 = zeros(1,4);
ky3 = zeros(1,4);
ky4 = zeros(1,4);

% RK4 coefficients
b = [1 2 2 1];

% Iterate, computing each K value in turn, then the i+1 step values
for i = 1:(N-1)        
    ky1(1) = dy1(t(i), y1(i), y2(i), y3(i), y4(i));
    ky2(1) = dy2(t(i), y1(i), y2(i), y3(i), y4(i));
    ky3(1) = dy3(t(i), y1(i), y2(i), y3(i), y4(i));
    ky4(1) = dy4(t(i), y1(i), y2(i), y3(i), y4(i));

    ky1(2) = dy1(t(i) + (h/2), y1(i) + (h/2)*ky1(1), y2(i) + (h/2)*ky2(1), y3(i) + (h/2)*ky3(1), y4(i) + (h/2)*ky4(1));
    ky2(2) = dy2(t(i) + (h/2), y1(i) + (h/2)*ky1(1), y2(i) + (h/2)*ky2(1), y3(i) + (h/2)*ky3(1), y4(i) + (h/2)*ky4(1));
    ky3(2) = dy3(t(i) + (h/2), y1(i) + (h/2)*ky1(1), y2(i) + (h/2)*ky2(1), y3(i) + (h/2)*ky3(1), y4(i) + (h/2)*ky4(1));
    ky4(2) = dy4(t(i) + (h/2), y1(i) + (h/2)*ky1(1), y2(i) + (h/2)*ky2(1), y3(i) + (h/2)*ky3(1), y4(i) + (h/2)*ky4(1));

    ky1(3) = dy1(t(i) + (h/2), y1(i) + (h/2)*ky1(2), y2(i) + (h/2)*ky2(2), y3(i) + (h/2)*ky3(2), y4(i) + (h/2)*ky4(2));
    ky2(3) = dy2(t(i) + (h/2), y1(i) + (h/2)*ky1(2), y2(i) + (h/2)*ky2(2), y3(i) + (h/2)*ky3(2), y4(i) + (h/2)*ky4(2));
    ky3(3) = dy3(t(i) + (h/2), y1(i) + (h/2)*ky1(2), y2(i) + (h/2)*ky2(2), y3(i) + (h/2)*ky3(2), y4(i) + (h/2)*ky4(2));
    ky4(3) = dy4(t(i) + (h/2), y1(i) + (h/2)*ky1(2), y2(i) + (h/2)*ky2(2), y3(i) + (h/2)*ky3(2), y4(i) + (h/2)*ky4(2));

    ky1(4) = dy1(t(i) + h, y1(i) + h*ky1(3), y2(i) + h*ky2(3), y3(i) + h*ky3(3), y4(i) + h*ky4(3));
    ky2(4) = dy2(t(i) + h, y1(i) + h*ky1(3), y2(i) + h*ky2(3), y3(i) + h*ky3(3), y4(i) + h*ky4(3));
    ky3(4) = dy3(t(i) + h, y1(i) + h*ky1(3), y2(i) + h*ky2(3), y3(i) + h*ky3(3), y4(i) + h*ky4(3));
    ky4(4) = dy4(t(i) + h, y1(i) + h*ky1(3), y2(i) + h*ky2(3), y3(i) + h*ky3(3), y4(i) + h*ky4(3));

    %y1(i+1)
    y1(i+1) = y1(i) + (h/6)*sum(b.*ky1);       
    y2(i+1) = y2(i) + (h/6)*sum(b.*ky2);       
    y3(i+1) = y3(i) + (h/6)*sum(b.*ky3);
    y4(i+1) = y4(i) + (h/6)*sum(b.*ky4);
end

txyz = [t,y1,y2,y3]

Answer 1

等待大神答复