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

% Function
    % Clarity z(1) = x; z(2) = y; z(3) = px; z(4) = py;
F = @(t, z) [z(3)/m; % X = q1
        z(4)/m; % Y = q2
        -2*k*z(1)+k*L*z(1)*(z(1)^2+z(2)^2)^(-.5)+k*z(1)*L*(z(1)^2+(z(2)-D)^2)^(-.5); % p1
        -2*k*z(2)+D*k+k*L*z(2)*(z(1)^2+z(2)^2)^(-.5)+k*L*(z(2)-D)*(z(1)^2+(z(2)-D)^2)^(-.5)-m*g]; % p2
% Jacobian Derivative Matrix: DF at y
Df = @(z) [0, 0, 1/m, 0;
    0, 0, 0, 1/m;
    -2*k+k*L*(z(1)^2+z(2)^2)^(-1/2)-k*z(1)^2*L*(z(1)^2+z(2)^2)^(-3/2)+k*L*(z(1)^2+(z(2)-D)^2)^(-1/2)-k*z(1)^2*L*(z(1)^2+(z(2)-D)^2)^(-3/2), -k*z(1)*L*z(2)*(z(1)^2+z(2)^2)^(-3/2) - k*z(1)*L*(z(2)-D)*(z(1)^2+(z(2)-D)^2)^(-3/2), 0, 0;
    -k*z(1)*L*z(2)*(z(1)^2+z(2)^2)^(-3/2) - k*z(1)*L*(z(2)-D)*(z(1)^2+(z(2)-D)^2)^(-3/2), -2*k+k*L*(z(1)^2+z(2)^2)^(-1/2)-k*z(2)^2*L*(z(1)^2+z(2)^2)^(-3/2)+k*L*(z(1)^2+(z(2)-D)^2)^(-1/2)-k*(z(2)-D)^2*L*(z(1)^2+(z(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,y) [F(t,y(1:4,:));
    Df(y(1:4,:))*y(5:8,:);
    Df(y(1:4,:))*y(9:12,:);
    Df(y(1:4,:))*y(13:16,:);
    Df(y(1:4,:))*y(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

