function stated = spring2Dn(t,state,n,m,k,b,x0,y0,xn,yn,Lr)

% (0,0)---(x1,y1)-...-(xm,ym)---(1,0)
% +x to the right, -y down
%
% state is [x; y; xd; yd]

% unpack state vector into x, y, xd, yd
x = state(1:n);
y = state(n+1:(2*n));
xd = state((2*n)+1:(3*n));
yd = state((3*n)+1:(4*n));

% add the anchors (x0,y0) and (xn,yn)
xx = [x0; x(:); xn];
yy = [y0; y(:); yn];

% compute lengths of each chain segment
L = zeros(n+1,1);
for i=1:n+1
   L(i) = sqrt((xx(i+1)-xx(i)).^2 + (yy(i+1)-yy(i)).^2); 
end
% compute stretch of each segment relative to resting length
S = (L-Lr).*((L-Lr)>0); % elastic bands; only pull, no push

g = 9.81;

% compute accelerations xdd and ydd of each state
xdd = zeros(n,1);
ydd = zeros(n,1);
for i=1:n
    cos1 = (xx(i+1)-xx(i))/L(i);
    cos2 = (xx(i+2)-xx(i+1))/L(i+1);
    sin1 = (yy(i+1)-yy(i))/L(i);
    sin2 = (yy(i+2)-yy(i+1))/L(i+1);
    xdd(i) = (-k*S(i)*cos1 + k*S(i+1)*cos2 - b*xd(i)      ) / m;
    ydd(i) = (-k*S(i)*sin1 + k*S(i+1)*sin2 - b*yd(i) - m*g) / m;
end

stated = [xd; yd; xdd; ydd];

end

