Clear["Global`*"]
ϵ=(1/4);
σ=1;
u[r_]:=(4*ϵ)*(((((σ/r)))^12)-((((σ/r)))^6)); (* Energy of two Argon atoms  as a function of the atom-atom distance. In dimensionless units of distance. *)

(* Question (a) *)
Plot[u[r],{r,0.7,1.5},AxesLabel->{"r","U(r)"},GridLines->Automatic,PlotRange->{-100, 100}] 
(*AxesLabel, GridLines, and PlotRange are absolutely neccesary;these make the reusltimg plot acceptable. Find good values for plotrange (placeholder values are used now) so that the plot looks similar to what's shown in the problem setup, with the key features, including the minimum,  clearly seen.  The example below clearly has the wrong range, needs to be changed.  *)

(*Let,(Subscript[r, 12]^( * )) = ???  You must find one that is appropriate, and explain why. Wrong placeholder value below. *)
r0= -100;

(* Question (b) *)

(*Using hard-coded Newton’s method*)
d1=D[u[y],{y,1}];
d2=D[u[y],{y,2}];
x0=r0;
y=x0;
x1=x0-d1/d2;
y=x1;
x2=x1-d1/d2;
y=x2;
x3=x2-d1/d2;
y=x3;
x4=x3-d1/d2;

Print["Starting Point, x0 = ",x0,"     F(x0)= ",u[x0]]
Print["First step,     x1 = ",x1,"     F(x1)= ",u[x1],"     Residual= ", x1-x0]
Print["Second step,    x2 = ",x2,"     F(x2)= ",u[x2],"     Residual= ", x2-x1]
Print["Third step,     x3=  ",x3,"     F(x3)= ",u[x3],"     Residual= ", x3-x2]
Print["Fourth step,    x4 = ",x4,"     F(x4)= ",u[x4],"     Residual= ", x4-x3]

(*Convergence Plot*)
ListLinePlot[{{0,u[x0]},{1,u[x1]},{2,u[x2]},{3,u[x3]},{4,u[x4]}}]

(* Question (c) *)
(*Using hard-coded Newton’s method*)
(* You need to run the random number generator it more than once to see the problem; by pure luck, some random starting points may be OK. It is also find to just stick in some randon values for x0 to see and explain what is going on.   *) 
x0=RandomReal[{-1.5,1.5}];
y=x0;
x1=x0-d1/d2;
y=x1;
x2=x1-d1/d2;
y=x2;
x3=x2-d1/d2;
y=x3;
x4=x3-d1/d2;

Print["Starting Point, x0 = ",x0,"     F(x0)= ",u[x0]]
Print["First step,     x1 = ",x1,"     F(x1)= ",u[x1],"     Residual= ", x1-x0]
Print["Second step,    x2 = ",x2,"     F(x2)= ",u[x2],"     Residual= ", x2-x1]
Print["Third step,     x3=  ",x3,"     F(x3)= ",u[x3],"     Residual= ", x3-x2]
Print["Fourth step,    x4 = ",x4,"     F(x4)= ",u[x4],"     Residual= ", x4-x3]

(*Convergence Plot*)
ListLinePlot[{{0,u[x0]},{1,u[x1]},{2,u[x2]},{3,u[x3]},{4,u[x4]}}]

Plot[u[r],{r,-1.5,1.5},AxesLabel->{"r","U(r)"},GridLines->Automatic,PlotRange->{-0.3,0.3}]

(* Question (d) *)
(* Old fashioned pencil and paper *)

(* After this point, for the solutions of the later problems, using NMinimize is the preferred solution. NMinimize attempts to find the global optimum. Other solutions using FindMinimum, and analytical solutions are acceptable, and applicable in some special cases; agreement between two different solutions is a good sanity check that the solution is correct.  Also analytical solution can be employed to check if your numerical solution is correct.*)
(*NMinimize is used here in 3 forms. The syntaxes are as follows:
NMinimize[f,x] minimizes f numerically with respect to x.
NMinimize[f,{x,y,…}] minimizes f numerically with respect to x,y,….
NMinimize[{f,cons},{x,y,…}] minimizes f numerically subject to the constraints cons.*)

(* Question (e) *)

(* Pencil and paper or simple code *)

(* Question (f) *)

;(* A general picture of the 3 atom molecule *)
(* With three atoms, follow the description of the problem; the energy function in question is (Subscript[U, tot])=(Subscript[U, 12])+(Subscript[U, 13])+(Subscript[U, 23]) *)
(Subscript[U, tot])[r1_,r2_,r3_]:=u[r1]+u[r2]+u[r3];

(* The numerical approach... *)
(* Numerically minimize (Subscript[U, tot]) for the three atom-atom distances... *)
NMinimize[{(Subscript[U, tot])[r1,r2,r3],r1>0&&r2>0&&r3>0},{r1,r2,r3}](* Using NMinimize, with constrains. Explain why we need these specific constraints r>0 here?  *)
FindMinimum[(Subscript[U, tot])[r1,r2,r3],{{r1, 100},{r2,100},{r3,100}}](* With FindMinimum, starting at (100,100,100) as a place holder starting point. That needs to be changed to a reasonable one. A good check if both NMinimize and FindMinimum give you the same answer.  *)

(* Note that FindMinimum is a _local_ minimizer, while NMinimize is a _global_ minimizer *)

(* Question (g) *)
(* If atoms are constrained to lie on a single line, we can redo the optimization in (f) with the associated constraints *)

;(* The key constraint, which we can identify from the image to the left, is r13=r12+r23 *)
(Subscript[U, totlin])[r1_,r2_]:=u[r1]+u[r2]+u[r1+r2]; (* We can rewrite the total energy function *)

(* Numerical Approach *)
NMinimize[(Subscript[U, totlin])[r1,r2],{r1,r2}](* Using NMinimize... *)
FindMinimum[(Subscript[U, totlin])[r1,r2],{{r1,100},{r2,100}}] (* Using FindMinimum. Starting point (100,100) is a wrong placeholder, you need to use a different, well justtified one *)

(* Extra Credit Question *)
(* Four atoms do not necessarily lie on a plane, so the native state is more complex than (f) *)
(* We must minimize the following energy function *)
U4[r12_,r13_,r14_,r23_,r24_,r34_] := u[r12]+u[r13]+u[r14]+u[r23]+u[r24]+u[r34];