function f = myforce_JMS(r,l,C,D)
% Force function for problem 2 of control 3, May.2022
% Input:
%   r(3,2) : particle positions (m)
%   l      : internal force parameter (m) 
%   C      : internal force parameter (N/m^3)
%   D      : external force parameter (N/m^2)
% Output:
%   f(3,2) : particle forces (N)

rab = r(:,2)-r(:,1);
fab = -C*(norm(rab)^2-l^2)*rab;
fext = -D*vecnorm(r).*r;
f = [-fab,fab] + fext;

end
