function [x, DFx] = lorenz_time1map(xo, s,b,r)
% evolve Lorenz flow for 1 time unit, including finding the jacobean matrix

F = @(t, y) [ -s*y(1,:)+s*y(2,:); r*y(1,:)-y(2,:)-y(1,:)*y(3,:); ...
              -b*y(3,:)+y(1,:)*y(2,:)];   % vec flow field
Df = @(y) [-s s 0; r-y(3,:) -1 -y(1,:); y(2,:) y(1,:) -b]; % DF at y

J0 = eye(3);   % initial Jac matrix

% 12-component ODE flow given by 3 components of solution and 9 components
% of the J matrix (J satisfies the ODE dJ/dt = Df.J)
G = @(t,z) [F(t,z(1:3,:)); Df(z(1:3,:))*z(4:6,:); Df(z(1:3,:))*z(7:9,:); ...
            Df(z(1:3,:))*z(10:12,:) ];

[ts, xs]  = ode45(G, [0 1], [xo; J0(:)],odeset('abstol',1e-3));    % numerically solve in t domain
x = reshape(xs(end,1:3), [3,1]);     % extract the answer at the final time t=1
J = xs(end,4:end);   % same for the J components
DFx = reshape(J,[3 3]);    % send J out as a 3x3 matrix.
