function [x, DFx] = Elastic_Pend_time1map(xo,g,L,k,m, D)
% Evolves Single elastic pendulum flow for 1 time unit
% Finds the jacobean matrix also.

% Function
F = @(t, y) [y(3)/m; 
    y(4)/m; 
    -k*y(1)+k*L*y(1)*(y(1)^2+(y(2)-D)^2)^(-1/2);
    -k*(y(2)-D)+k*L*(y(2)-D)*(y(1)^2+(y(2)-D)^2)^(-1/2)- m*g];
% Jacobian Derivative Matrix: DF at y
Df = @(y) [0, 0, 1/m, 0;
    0, 0, 0, 1/m;
-k + k*L*(y(1)^2+(y(2)-D)^2)^(-1/2)-k*L*y(1)^2*(y(1)^2+(y(2)-D)^2)^(-3/2), -k*L*y(1)*(y(2)-D)*(y(1)^2+(y(2)-D)^2)^(-3/2), 0, 0;
-k*L*y(1)*(y(2)-D)*(y(1)^2+(y(2)-D)^2)^(-3/2), -k+k*L*(y(1)^2+(y(2)-D)^2)^(-1/2)-(y(2)-D)^2*k*L*(y(1)^2+(y(2)-D)^2)^(-3/2), 0, 0]; 

J0 = eye(4);   % Initial Jacobians matrix

% 16-component ODE flow: 4 components of solution and 12 components for 
% the J matrix. (J satisfies the ODE dJ/dt = Df.J)
G = @(t,z) [F(t,z(1:4,:));
    Df(z(1:4,:))*z(5:8,:);
    Df(z(1:4,:))*z(9:12,:);
    Df(z(1:4,:))*z(13:16,:);
    Df(z(1:4,:))*z(17:20,:) ];

% Numerically solves in time domain
[ts, xs]  = ode45(G, [0 1], [xo; J0(:)],odeset('abstol',1e-3)); 
% Extracts the solution at the final time t=1.
x = reshape(xs(end,1:4), [4,1]); 
J = xs(end,5:end);   % Extracts J components
DFx = reshape(J,[4 4]);    % send J out as a 4x4 matrix.

end

