function E = Potencial_C1_JMS(d,d0,V0)
% Solution to problem 2 of Control 2.
% Input:
%   d  : distance(s) (nm)
%   d0 : distance of the energy minimum (nm)
%   V0 : depth of the energy minimum (eV)
% Output:
%   E  : potential energy for each distance (eV)
% JMS, Oct.2022

% Internal parameters
d1 = d0*(1+sqrt(3));
a = V0/d0^2;
b = -2*a/3/(d1-d0);

% Energy
E = zeros(size(d));
n = numel(d);
for i = 1:n
  if d(i)<=d0
    E(i) = a*(d(i)^2-2*d0*d(i));
  elseif d0<d(i) && d(i)<=d1
    E(i) = a*(d(i)^2-2*d0*d(i)) + b*(d(i)-d0)^3;
  end
end

% A faster way
% E = a*(d.^2-2*d0*d) .* (d<=d1) + b*(d-d0).^3 .* (d0<d & d<=d1);

end
