Correcting matrix dimensions in Runge-Kutta ODE solver

I have a simple SIR model written in a .m file:

function dx = test1(t, y)
dx = [0; 0; 0];
beta = .003; 
delta = 1;
dx(1) = -beta * y(1) * y(2);
dx(2) = beta * y(1) * y(2) - delta * y(2);
dx(3) = delta * y(2);
end

And I made a Runge-Kutta ODE solver:

function [t y] = ode_RK4(f,dt,t_end,y0,steps)

h = (t_end - dt) / steps;
h_2 = h / 2;
y(1,:) = y0;
t(1) = dt;
h_6 = h/6;

for i = 1 : steps
  t(i+1) = t(i) + h;
  th2 = t(i) + h_2;
  k1 = f(t(i), y(i,:));
  k2 = f(th2, y(i,:) + h_2 * k1);
  k3 = f(th2, y(i,:) + h_2 * k2);
  k4 = f(t(i+1), y(i,:) + h * k3);
  y(i+1,:) = y(i,:) + (k1 + k2+k2 + k3+k3 + k4) * h_6;
end
end

Now when I want to solve the SIR model using

 [t,y] = ode_RK4(@test1,0,10,[1000 0 1],10)

I get an Error using +. Matrix dimensions must agree. Error in ode_RK4 (line 43). I am quite not sure where to make corrections in the RK solver and would be very grateful if anyone can help me to solve my problem.

[Torsten]

It might suffice do remove the line

dx = [0; 0; 0];

in "test1".

Best wishes

Torsten.

[Denis Benka]

I tried to remove the dx line in test1.m but then I got this output from the solver

 y =

        1000           0           1
        1000           0           1
        1000           0           1
        1000           0           1
        1000           0           1
        1000           0           1
        1000           0           1
        1000           0           1
        1000           0           1
        1000           0           1
        1000           0           1

[Torsten]

Since the initial condition for y(2) is y(2)=0, this is the correct result.

Best wishes

Torsten.

[Denis Benka]

I deeply apologize, my bad... I did not realize that, thank you! Now when I put the initial condition in correct order [1000 1 0] it works. Thank you again and have a nice day!