function f=rhslyap(x,y)
global b r sigma
n=3;
%   Integrate phase point (y(1)=x, y(2)=y, y(3)=z)
dydx(1)=sigma*(y(2)-y(1));
dydx(2)=r*y(1)-y(1)*y(3)-y(2);
dydx(3)=y(1)*y(2)-b*y(3);
%   Integrate tangent vectors with Jacobian (first component myin, second
%   component myin+1, third component myin+2, for myin =4,7,10);
%   for Lorenz the Jacobian is   |-sigma   sigma    0|
%                                | (r-z)    -1     -x|
%                                |   y       x     -b|
myin=n+1:n:n*(n+1);
dydx(myin)=-sigma*y(myin)+sigma*y(myin+1);
dydx(myin+1)=(r-y(3))*y(myin)-y(myin+1)-y(1)*y(myin+2);
dydx(myin+2)=y(2)*y(myin)+y(1)*y(myin+1)-b*y(myin+2);
f=dydx';