function rhs = lorentz_rhs(t, y, dummy, sigma, rho, beta)
rhs = [sigma*(y(2) - y(1))
       y(1)*(rho - y(3)) - y(2)
       y(1)*y(2) - beta*y(3)];