function yprime = lorenzode(t,y,b,s,r)
% Lorenz differential equations
yprime = [s*(y(2)-y(1)) ; r*y(1) - y(2) - y(1)*y(3); y(1)*y(2) - b*y(3)];