Unable to let Ringe Kutta 4 two ODE work

Hi all,

I do not understand why my ringe kutta 4 method won't work. Maybe it is a small fault.

the differential equations are:

dy1/dx = (a - alpha*y2)*y1 - W(t)*y1,

dy2/dx (-c+gamma*y1)*y2

Can please somebody help me???

clear all
close all
clc

% parameter
a = 3;
c = 3;
alpha = 2*10^(-3);
gamma = 7*10^(-3);

% define time
ti = 0;
tf = 5;
delta_t0 = 1;

t0 = ti:delta_t0:tf;
h0 = delta_t0;
N0 = ceil(tf/h0);

% define the seasonal fishing load factor
W = zeros(N0,1);

for n =1:N0
W(n,1)=fishing_load_factor(t0(n));
end


% Define functions
fy1 = @(t,y1,y2,W) (a-alpha*y2)*y1-W*y1;
fy2 = @(t,y1,y2) (-c+gamma*y1)*y2;

%initial conditions
y1(1)=600;
y2(1)=1000;
t0(1)=0;

% Update loop
for i = 1:N0
    %Update time
    t0(1+i)=t0(i)+h0;
    
    % Update y1 and y2
    K11 = fy1(t0(i)      ,y1(i)          , y2(i), fishing_load_factor(t0(i)));
    K12 = fy2(t0(i)      ,y1(i)          , y2(i));
        
    K21 = fy1(t0(i)+h0/2  ,y1(i)+h0/2*K11, y2(i)+h0/2*K12,  fishing_load_factor(t0(i)+h0/2));
    K22 = fy2(t0(i)+h0/2  ,y1(i)+h0/2*K11, y2(i)+h0/2*K12);
    
    K31 = fy1(t0(i)+h0/2  ,y1(i)+h0/2*K21, y2(i)+h0/2*K22, fishing_load_factor(t0(i)+h0/2));
    K32 = fy2(t0(i)+h0/2  ,y1(i)+h0/2*K21, y2(i)+h0/2*K22);
 
    K41 = fy1(t0(i)+h0    ,y1(i)+h0*K31  , y2(i)+h0*K32, fishing_load_factor(t0(i)+h0));
    K42 = fy2(t0(i)+h0    ,y1(i)+h0*K31  , y2(i)+h0*K32);
    
    y1(1+i) = y1(i)+h0/6*(K11 + 2*K21 + 2*K31 + K41);
    y2(1+i) = y2(i)+h0/6*(K12 + 2*K22 + 2*K32 + K42);
    
    
end