The first problem is that you begin by assigning to mu after which you define the equations. Thereby the equations won't change if you change mu later.
The next problem is that even if you postpone assigning to mu till you use dsove/numeric then that sol will be the solution corresponding to the assigned value of mu.

The solution is to use the parameters option in dsolve/numeric, see
?dsolve,numeric,interactive

So with the necessary changes you got:

restart;
with(plots):
with(DEtools):
k := 100; L := .25; M := .5; g := 4;
F[N] := k*L*(1-L/sqrt(x(t)^2+M^2));
eq1 := diff(x(t), t) = v(t);
eq2 := diff(v(t), t) = k*x(t)*(1-L/sqrt(x(t)^2+M^2))+mu*abs(F[N]-g)*signum(x(t));
sol := dsolve({eq1, eq2, v(0) = 1, x(0) = .1098},numeric,parameters=[mu]);
sol(parameters=[0.2e-1]); #Here mu will be replaced by that value
newplots := [odeplot(sol, [t, x(t)], t = 0 .. 2, colour = red), odeplot(sol, [x(t), v(t)], t = 0 .. 2, colour = blue)];
display(array(newplots));
#Now your loop:
for i to 17 do
    mu := (i-9)*(1/2); #now assignment is harmless
    sol(parameters=[mu]);
    sss := cat("The value of μ is ", convert(mu, string), ".");
    P[i] := odeplot(sol, [t, x(t)], t = 0 .. 2, colour = red, title = sss)
end do:

display(seq(P[k], k = 1 .. 17), insequence = true);

Using equation 3.4.82 in the reference you gave and also
diff(phi,t,t) = phidot*diff(phidot,phi) we are back to the second order ode:

ode:=diff(phi(t),t,t)=-(sqrt(3/2)/M^2*sqrt(diff(phi(t),t)^2+m^2*phi(t)^2)*diff(phi(t),t)+m^2*phi(t));
#Following Adri van der Meer we turn the equation into a first order system:
sys:=diff(phi(t),t)=phi1(t),subs(diff(phi(t),t)=phi1(t),ode);
#I pick my own parameters, m = M = 1:
with(DEtools):
DEplot(eval([sys],{M=1,m=1}),[phi(t),phi1(t)],t=0..30,[[phi(0)=1,phi1(0)=0],[phi(0)=-1,phi1(0)=0],[phi(0)=2,phi1(0)=-1],[phi(0)=-2,phi1(0)=1]],color=grey,linecolor=blue,stepsize=.1,thickness=1);

You got a sequence of 1x1 matrices!
Assign that sequence to L.
Then
V:=Vector([L]);
will be a column vector with the elements from the 1x1 matrices.
Now export with ExportVector or with ExcelTools:-Export.

I don't quite see in your code any attempt to implement the piecewise defined A at the top. Presumably (t,365) is meant as t-floor(t/365)*365.

Knowing nothing about stochastic differential equations the following should be taken as an experiment.
The example used is only somewhat similar to yours.

Added: A somewhat simpler looking way at the bottom.

restart;
#The ode considered is
ode:=diff(U(t),t)=-(Ap(t)+r(t)+B*U(t))*U(t);
#where Ap is a piecewise defined function, B is a constant, and r is a random variable here defined as:

r:=RandomTools:-Generate(distribution(Normal(0, .5)), makeproc=true);
#We shall not use ode at all, but use procedural input to dsolve/numeric:

p:=proc(N,t,Y,YP) local Ap,B;
   B:=.1;
   Ap:=piecewise(t<=1,.1+ 0.4*t,0.5+ 0.004*t); 
   YP[1]:=-(Ap+r(t)+B*Y[1])*Y[1]
end proc;
#It appears to me that the only reasonable method is a forward method with no error corrections.
#Thus I use the simple Euler method classical[foreuler].
#Choose some stepsize consistent with the application we have in mind.
#Here chosen arbitrarily to 0.01.
res:=dsolve(numeric,procedure=p,number=1,initial=Array([1]),start=0,procvars=[U(t)],method=classical[foreuler],stepsize=0.01);
plots:-odeplot(res,[t,U(t)],0..2);

#########################################
##Another way and I guess simpler looking.
restart;
ode:=diff(U(t),t)=-(Ap(t)+r(t)+B*U(t))*U(t);
R:=RandomTools:-Generate(distribution(Normal(0,.5)), makeproc=true);
r:=proc(t)
   if not type(t,numeric) then return 'procname(_passed)' end if;
   R()
end proc;
Ap:=t->piecewise(t<=1,.1+ 0.4*t,0.5+ 0.004*t);
B:=0.1:
res:=dsolve({ode,U(0)=1},numeric,known=r(t),method=classical[foreuler],stepsize=0.01);
plots:-odeplot(res,[t,U(t)],0..2);

Taking jj=2 and keeping the rather large g:

soln:=dsolve({initcondNORMAL, op(subs(k = 1, n = 2, g = 10, [eq1,eq2,eq3]))}, numeric,output=listprocedure);
Z,Z1,X,X1:=op(subs(soln,[z(t),diff(z(t),t),x(t),diff(x(t),t)]));
plots:-odeplot(soln,[t,z(t)],0..0.15);
fsolve(Z1); #The time for maximal z
Z(%); #The value of zmax
T:=fsolve(Z,.15); #The time for z = 0.
arctan(Z1(T)/X1(T)); #The angle in radians
%*180/evalf(Pi); #The angle in degrees
####################
Try also

solE:=dsolve({initcondNORMAL, op(subs(k = 1, n = 2, g = 10, [eq1,eq2,eq3]))}, numeric,events=[[diff(z(t),t),halt],[z(t),halt]]);
solE(.1);
solE(eventclear);
solE(.2);

If you can live with a parenthesis aound the exponent then beta=10^``(-6) is a solution.

A rather weird way:

ten:=convert(10,symbol);
ten^(-6); #Doesn't do it, but:
series(ten^(-6),ten);

Although your example is simple, it is a polynomial, so fsolve returns 2 solutions (if they are real).
So I have taken the first one found (the [1] attached to fsolve):
p:=(a,b)->fsolve( x^2 + a*x + b =0, x)[1];
#Use argument free syntax:
plot3d( p , 0.1..5, 0.5..3,axes=boxed,labels=[a,b,p]);

When you don't have a polynomial fsolve only returns one solution, so omit the [1].
Example:

p:=(a,b)->fsolve( x + a*ln(x) + b =0, x);
plot3d( p , 0.1..5, 0.5..3,axes=boxed,labels=[a,b,p]);

You may want to take a look at the events option:

?dsolve,events

restart;
sys:=diff(x(t),t)= x(t)+y(t) , diff(y(t),t)= x(t)-y(t);  
res:=dsolve({sys,x(0)=1/2,y(0)=0},numeric,events=[[x(t)-1,halt],[y(t)-2,halt]]);
plots:-odeplot(res,[x(t),y(t)],0..5);
plots:-odeplot(res,[t,x(t)],0..5);

Here is a simple example, which probably just shows that I don't understand your problem:

restart;
with(LinearAlgebra):
M:=<seq(m[i](t),i=1..3)>;
H:=<M[2],M[1],M[3]>; #Simple example
eqs:=diff~(M,t)+ M &x H-a*M &x diff~(M,t)/DotProduct(M,M,conjugate=false)=~ <0,0,0>;
sys:={seq(eqs[i],i=1..3)};
res:=dsolve(eval(sys,a=1) union {m[1](0)=1,m[2](0)=0,m[3](0)=0}, numeric);
plots:-odeplot(res,[seq([t,m[i](t)],i=1..3)],0..10);

Here is a simple example:

v:=Vector(datatype=float);
for i from 1 to 10 do v(i):=subs(Optimization:-NLPSolve(sin(x)*x^i, x=1..5,method=branchandbound,maximize)[2], x) end do:
v;

Then use ExportVector.

In your infinity version there appear factors beta(beta), they should just be beta. See comment at the end, though.
Specifically, you have

1/c1 = 6.274557048*10^8*(Int((32.*sin(0.5000000000e-1*beta)/beta+(960.*(cos(0.5000000000e-1*beta)-40.*sin(0.5000000000e-1*beta)/beta+1600.*sin(0.2500000000e-1*beta)^2/beta^2))/beta^2)^2/((75*(1.50*tanh(beta)*beta(beta)+1))/(75+0.2e-1*tanh(beta)*beta(beta))+2.32*coth(beta)*beta(beta)), beta = 0 .. infinity));
#You see the appearance of beta(beta).
#Let A be the integral (using Int):

A:=Int((32.*sin(0.5000000000e-1*beta)/beta+(960.*(cos(0.5000000000e-1*beta)-40.*sin(0.5000000000e-1*beta)/beta+1600.*sin(0.2500000000e-1*beta)^2/beta^2))/beta^2)^2/((75*(1.50*tanh(beta)*beta(beta)+1))/(75+0.2e-1*tanh(beta)*beta(beta))+2.32*coth(beta)*beta(beta)), beta = 0 .. infinity);
#Now do
subs(beta(beta)=beta,A);
B:=normal(%);
#This integral can be evaluated, but it is not all that easy. The problem is not at beta = 0.
f:=IntegrationTools:-GetIntegrand(B);
asympt(f,beta, 9);
series(f,beta=0, 9);
plot(f,beta=0..100);
evalf(Int(f,beta=0..1000));
evalf(Int(f,beta=0..5000));

More work could be done on this, of course. You could try something like:

evalf[15](Int(f,beta=0..10000,method = _d01akc));
#and
evalf[15](Int(f,beta=0..200000,method = _d01akc,epsilon=1e-10));
#for large values of beta f is bounded from above by
g:=eval(f,{sin=1,cos=1,tanh=1,coth=1,-38400.=38400});
#so you can get an upper bound for the error by computing:
evalf[15](Int(g,beta=200000..infinity));

Comment on beta(beta):

Try the following
restart;
evalf(Int(beta(beta),beta=0..1));
evalf(Int(beta,beta=0..1));
int(beta(beta),beta=0..1);
#The reason the numerical integration works is that in Maple the result of
5(5);
#is 5, as is
5(x);

restart;
g:=y*(1-1/(x^2+y^2)):
# g = 0.1
plots:-contourplot(g,x=-1..1,y=0..2,contours=[0.1]); p:=%:
dt:=plottools:-getdata(p):
L:=ListTools:-Flatten(select~(type,[dt],Matrix)):
M:=<op(L)>;
#plot(M);
#M[1..10,..];
#Lots of repeats, so:
Mu:=Matrix(ListTools:-MakeUnique(convert(M,listlist)));
#plot(Mu);

Then do the same for g = 0 and save in Mu0, say. Notice now that plot(Mu0); connects points in a confusing manner. But you can do plot(Mu0,style=points) to avoid seeing these connecting lines.
To export Mu and Mu0 to Excel see ?ExcelTools

phi[i], i=1..N are not defined as functions (procedures) so should not be used as such. 

restart;
N:=4:
for i from 1 to N do
 phi[i]:=piecewise(t>=x[i-1] and x[i]>t, (t-x[i])/(h), x[i]<t and t<=x[i+1], (x[i+1]-t)/(h), 0);
end do;
V:=Int(diff(phi[1], t)*diff(phi[1], t), t=x[0]..x[2]);
value(V) assuming x[0]<x[1],x[1]<x[2];

Since log(y) - > -infinity as y - > 0+, there wouldn't be room on your screen for your plot. So a compromise is made.

Try for comparison

plots:-loglogplot(f, x = 0 .. 5*10^10, y = 1/1000 .. 2, thickness = 3, discont = true);

restart;
#Taking Alejandro Jakubi's example:
S:= A*cos(2*x)*sin(3*z) + B*cos(7*x)*sin(4*z) + C*sin(3*z)+K*sin(4*z)+cos(7*x);

#Using frontend:
eval(frontend(coeff,[S,sin(3*z)]),{cos=0,sin=0});
eval(frontend(coeff,[S,cos(7*x)]),{cos=0,sin=0});
########
Edited: I have kept my first version of a procedure for this task below, but here is a much shorter one using coeffs and no freezing:

p:=proc(S::`+`) local cs,T;
   cs:=[coeffs(S,indets(S,specfunc(anything,{sin,cos})),'T')];
   select(type,`[]`~([T],cs),[specfunc(anything,{sin,cos}),anything])
end proc;
p(S);

############
Original version:

Making a procedure to find all. This is in parts like my answer to
http://www.mapleprimes.com/questions/142564-Extracting-Coefficients-From-A-Fourier-Series

p:=proc(S::`+`) local L,Lf,cs,ff,ffcs;
   L:=convert(S,list);
   Lf:=evalindets(L,specfunc(anything,{sin,cos}),freeze);
   cs:=eval(L,{sin=1,cos=1});
   ff:=Lf/~cs;
   ffcs:=`[]`~(ff,cs);
   thaw(remove(type,ffcs,[`*`,anything]))
end proc;
   
p(S);