function f=lorenzrhs(t,y)
% RHS for Lorenz equations
global sigma r b
f=[sigma*(y(2)-y(1));r*y(1)-y(1)*y(3)-y(2);y(1)*y(2)-b*y(3)];