1. In procedure generate_ic, p2_sol is a sequence (at least with my Maple 2015 version) of two values containing p1.
   Inset this print line to see that:

   p2_sol := solve(H_eq, p2); 
   print('p2_sol'=p2_sol);
   

2. If p2_sol is indeed a sequence forn your Maple version, then this (simpler) command fails 

   select(x -> 0 < evalf(x), p2_sol);
   Error, invalid input: select expects 2 or more arguments, but received 1
   

   and you must write instead

   select(x -> 0 < evalf(x), [p2_sol]);  # note the square brackets
   

3. As p2_sol contains the name p1, no testing is applicable because Maple does not know anything about p1.
   To understand that run this simplified version of generate_ic:

   test_1 := proc (q1_val, ET_val, f_val) 
     local term, H_eq, p2_sol, p2_val; 
     term := sqrt(q1_val^2+1)-1+f_val; 
     H_eq := (1/2)*p1^2+(1/2)*p2^2-f_val+(1/2)*term^2 = ET_val; 
     p2_sol := [solve(H_eq, p2)]; 
     print('p2_sol'=p2_sol);
     select(x -> 0 < evalf(x), p2_sol)
   end proc:
   
   interface(displayprecision=5):
   
   test_1(.1, ET, f)
   
            p2_sol = [2.000*10^(-8)*sqrt(-2.500*10^15*p1^2+8.745*10^13), -2.000*10^(-8)*sqrt(-2.500*10^15*p1^2+8.745*10^13)]
   
   Error, (in test_1) selecting function must return true or false
   

Last but not least: can you provide a picture of your spring pendulum (is it a simple or double one or a mixed one?)
 

The worksheet is written with Maple 2015 and the output of solve(eqn1, allsolutions); (see @dharr ) differ from what you get in newer version.
This is why I explicitely split this output in two components named ans_2015_0 and ans_2015_1 wich corresponds to @dharr's ans1 and ans2 resoectively.

The plot serves as a guidance but is not at all mandatory: what matters is what comes next.

sets_mmcdara.mw

restart
solve(sqrt(2)*sin(2*x-(1/6)*Pi) = 1)
                             5    
                             -- Pi
                             24   
solve(sqrt(2)*sin(2*x-(1/6)*Pi) = 1, allsolutions=true)
                   5       3                
                   -- Pi - - Pi _B1 + Pi _Z1
                   24      4                
about(_B1)
Originally _B1, renamed _B1~:
  is assumed to be: OrProp(0,1)

about(_Z1)
Originally _Z1, renamed _Z1~:
  is assumed to be: integer

# BTW: you missed the Pi constant in your solution

eval(sqrt(2)*sin(2*x-(1/6)*Pi), x=5/24) = 1
                    (1/2)    /  5    1   \    
                  -2      sin|- -- + - Pi| = 1
                             \  12   6   /    
eval(sqrt(2)*sin(2*x-(1/6)*Pi), x=5/24*Pi) = 1
                             1 = 1

Download Allsolutions.mw

Simply use  package Student:-Statistics  (it's a shame that package Statistics forgot ExploreRV !)
 

with( Student:-Statistics ):
X := NormalRandomVariable(mu, sigma);
ExploreRV( X );


# Or, even more simply:
Student:-Statistics:-ExploreRV()

Raw result

By the time I was writting this reply Scot Gould was submitting his one.
Nevertheless I decided not to delete mine, even if there is likely nothing new compared to Scot's.

Customize_it.mw

You can animate wrt 1 or more parameters by selecting them.

BTW: You are Salim Barzani aren't you? 

I guess there is a conflct netween background in G1 and background in animate.
Here is a Maple 2015 Workaround.mw

Feel free to adjust the plotting ranges in G1.

Here is another alternative (more funny?)

f := unapply(PDF(Normal(mu, sigma), x), x):

G3 := proc(t)
  plots:-display(
    DensityPlot(X, color = blue, thickness = 5, range = mu-4*sigma-1..t)
    , plot(f(tau), tau = mu-4*sigma .. t, thickness = 5, color = red)
    , gridlines=true
    , view=[mu-4*sigma .. t, default]
  )
end proc:

plots:-animate(G3, [t], t = mu-4*sigma .. mu+4*sigma, background="AliceBlue") 

Download stop-working_mmcdara.mw

Point 1:
The signature of [2, 1, 4, 3] (your example p2) is +1 (even permutation) and not -1!

Point 2:
You can indeed ask Geminy for a code, but you can also ask Mapleprimers, or more simply search for parity in the help pages.
This is what an AIA (Astute IA, see @nm's reply) should answer you

Download PermParity.mw

@nm: It would be interesting to see what the IA you used in your reply answers you for this parity problem: does it know about GroupTheory package?

I suggest you to read the help pages to see what the correct syntaxes are

Download LA.mw

The soution is based upon two procedures f and `w'/w` :
 

Download series-finding_mmcdara.mw

Download collection_of_coefficients_mmcdara.mw

An ordinary differential equation of order n writes 

diff(phi(x), x$n) = F(x, diff(phi(x), x), ..., diff(phi(x), x$(n-1))

Your differential equation (DE) is 

(diff(phi(x), x)^2 = -v(phi(x))

and is not an ODE.

Undortunately dsolve/numeric can handle only ODEs, not DEs.

 

One idea could be to rewrite your DE this way where siga is either -1 or +1 dependeing on the branch you want to select.

diff(phi(x), x = sigma * sqrt(-v(phi(x)))

Nevertheless this will not give you a correct result as this simple example proves:

   ( initial worksheet where @dharr detected a typo  Download DE_not_ODE.mw)

Edited worksheet DE_not_ODE_edited.mw

Every time you get an error while solving numerically a BVP yiu must look the dsolve[numeric][BVP] help page.
The Error, (in dsolve/numeric/bvp) initial Newton iteration is not converging is a classic for which this help page proposes a few cures (search the word Newton).
Another important help page is dsolve[numeric_bvp,advanced].

Here is the solution to yoiur problem (I just treated b1 and b2: end the job by adjusting the continuation term, maxmesh and abserr as you want.

VC1_mmcdara.mw