function [FspringNew, ST] = Sub1(N, Fspring, S, DeltaD, Fmax, Fmin)

FspringNew = zeros(N,1);
layerDeltaD = zeros(N,1);
% Set the tangent stiffness based on the original stiffness
ST = S;
for n = 1:N

    if n == N
        % This is the bottom layer so the displacement is the nodal
        % displacement
        layerDeltaD(n) = DeltaD(n);
    else
        % Otherwise we need to find the difference between the nodal
        % displacements
        layerDeltaD(n) = (DeltaD(n) - DeltaD(n+1));
    end
    FspringNew(n) = Fspring(n) + S(n)*layerDeltaD(n);
    %Check for yeilding
    if FspringNew(n) > Fmax(n)
        % We yielded so update the spring force to the maximum spring force
        FspringNew(n) = Fmax(n);
        % Calculate a stiffness that is compatible with the new spring
        % force
        ST(n) = (FspringNew(n)-Fspring(n))/layerDeltaD(n);
    elseif FspringNew(n) < Fmin(n)
        % We yielded so update the spring force to the minimum spring force
        FspringNew(n) = Fmin(n);
        % Calculate a stiffness that is compatible with the new spring
        % force
        ST(n) = (FspringNew(n)-Fspring(n))/layerDeltaD(n);
    end
 
end