z is a global variable in the first and second session. The assignment b:=2*z is saved.
In the third session z is declared local and assigned to k. The global version of b is read. w is assigned to b, which was assigned to the global version of z. Thus w will not be affected by assignments to the local z.

Try changing the declarations to
local w; global z;

The difference is exhibited here:

restart;
ilog[3](3^22-1); #output 22
restart;
Digits:=50:
ilog[3](3^22-1); #output 21

Try

showstat(ilog);

You will see that in the indexed case, using ilog[b], where the index b is 10 or 2, then the builtin functions ilog10 or ilog2 are used, respectively.
In all other indexed cases we have b1 := evalf(b); immediately introducing results possibly depending on the setting of Digits.

V:=1-2*M*exp(-2*l*r)+q^2*exp(-4*l*r);
plot(eval(V,{M=.1,l=.4,q=.7}),r=0..3,filled,color=blue);

Change
part:= int(list[i]*f*(x+1),x=-1..1)*list[i]:
into
part:= int(list[i](x)*f*(x+1),x=-1..1)*list[i](x):

since 'list' is a list of procedures (functions).

Then

F:=functions(3);
returns 3 functions and you can do e.g. F[3](5);

Alternatively (and better, I think) wait with unapply until the very end:

functions:= proc(n)
local L, list, p, f, sum, i, part, g, normg, x:
L:=1/sqrt(2):
list:=[L]:

for p from 2 to n do
  f := x**(p-1):
  sum := 0:

  for i from 1 to (p-1) do
    part:= int(list[i]*f*(x+1),x=-1..1)*list[i]:
    sum := sum + part:
  end do:

  g:= f-sum:
  normg:= sqrt(int(g^2*(x+1),x = -1..1)):
  L:=g/normg:
  list := [op(list), L]:
end do:
unapply~(list,x)
end proc:

You can now do as before.

Let u be your expression. Then try:

evalindets(u,radical,factor);
simplify(%) assuming positive;
q:=denom(%)/lambda;
simplify(algsubs(q =w,%%));
expand(%);
eval(%,w=q);
map(simplify,%) ;
combine(%) assuming positive;

I'm not following you. I cannot see that f1 and f2 are solutions:

odetest(y(t)=f(t),ode2); #Yes, y(t)=f(t) is a solution
odetest(y(t)=f(t-t1),ode2);
eval(%,{t=0,t1=1,s=1});
evalf(%);
#Not zero meaning that y(t)=f(t-t1) is not a solution of ode2.

However, reduction of order gives you a solution linearly independent of f(t) easily:

simplify(eval(ode2,y(t)=f(t)*v(t)));
eval(%,diff(v(t),t)=w(t));
dsolve(%);
combine(int(eval(rhs(%),_C1=1),t)); #Finding v(t)
f(t)*%;
convert(%,trigh);
f2:=simplify(%); #The second solution
odetest(y(t)=f2,ode2); #Indeed
#There is really no reason to check the wronskian since v(t) is not constant in t.
#But if we do:
VectorCalculus:-Wronskian([f(t),f2],t,determinant)[2];
convert(%,exp);
#we find (1/64)*s^2.
# s= 0 is a special case and not very interesting:
eval(ode2,s=0);

You cannot have infinity as the right "endpoint". Choose a large number in this case I chose 10.
I assume that your last boundary condition is meant to mean f''(infinity) = 0. With infinity = N that would be
(D@@2)(f)(N) =0 or (D@D)(f)(N)=0.
Try this then.
N := 10;
bcs := (D(f))(0) = 1, f(0) = 0, (D(f))(N) = 0, ((D@@2)(f))(N) = 0;
sol := dsolve({bcs, eq1}, numeric, method = bvp[midrich], maxmesh = 1024);
odeplot(sol,[eta,f(eta)],0..N);

u:=BesselJ(1,rho*exp(I*Pi/4));
diff(Re(u),rho); # Notice abs(1, ...)
eval(%,rho=0.5);
evalf(%);
# abs(1,w) is the derivative of abs(w)
#try
diff(abs(w(rho)),rho);
# What you are after may be what you get by postponing using Re:
u1:=diff(u,rho);
eval(u1,rho=0.5);
evalf(%);
Re(%);

There is probably no closed form available in the general case.
In two very special cases, yes:

dsolve(eval(ODE[4],a=0));
dsolve(eval(ODE[4],epsilon[2]=0));

Otherwise, solve numerically by maybe using the parameters option in dsolve/numeric as illustrated here.
You need initial conditions. I have used y = 0, but included the values as parameters Y0 and Y1:

res:=dsolve({ODE[4],Y(0)=Y0,D(Y)(0)=Y1},numeric,parameters=[Y0,Y1,a, beta, omega, epsilon[1], epsilon[2], mu[0]]);
#Now set values for the parameters. With this many parameters it is difficult to remember the order, but instead of giving numbers in the correct order you can give equations (in any order) as I have done here:

res(parameters=[Y0=1,Y1=0,a=.1, beta=.2, omega=.3, epsilon[1]=.4, epsilon[2]=.5, mu[0]=-6]);
#This would have been just as good:
res(parameters=[a=.1, beta=.2, omega=.3, epsilon[1]=.4, epsilon[2]=.5, mu[0]=-6,Y0=1,Y1=0]);
#Now plotting
plots:-odeplot(res,[y,Y(y)],0..25);

Q:=Matrix([[a+b,c+d],[e+f,g+h]]);
M:=Matrix([[A,B],[C,F]]);
assign~(M=~Q);
M;
A,B,C,F;

restart;
ode0:=diff(f(x),x,x)/((1+diff(f(x),x)^2)^(3/2)) = e*(x*sin(alpha)+f(x)*cos(alpha))+P;
ics:=D(f)(0)=0, f(0)=1;
params:={e=1/10,alpha=Pi/3};
ode:=eval(ode0,params);
res:=dsolve({ode,ics},numeric,output=listprocedure,parameters=[P]);
F,F1:=op(subs(res,[f(x),diff(f(x),x)]));
K:=1/sqrt(3); #So f'(a) = -1/sqrt(3)
#The procedure p surves dual purposes.
#The default is the primary: Finding f(a) and f'(a)+K.
#The secondary purpose is for use in exploring (plotting).
p:=proc(a,P,output::{procedure,list(procedure)}:='primary') local fa,f1a,i; global p1,p2;
   res(parameters=[P]);
   if output<>'primary' then return output end if;
   p1(_passed):=F(a);
   p2(_passed):=F1(a)+evalf(K);
   if type(procname,indexed) then cat(op(0,procname),op(procname))(_passed) else [p1,p2](_passed) end if
end proc;
p1:=proc(a,P) p[1](_passed) end proc;
p2:=proc(a,P) p[2](_passed) end proc;
#Hinting at exploration
plot(p(.9,-.7,[F,F1]),0..1.5);
aa:=1.9:
plots:-animate(plot,['p(aa,P,[F,F1])',0..aa,-K..1],P=-.7..-.06);
#Solving
sol:=fsolve([p1,p2],[1.9,-.6]);
plots:-odeplot(res,[[x,f(x)],[x,diff(f(x),x)]],0..sol[1]);
F(sol[1]);
F1(sol[1])+evalf(K);

I suggest that you give us the actual ode. In the meantime here is a very simple example, where it is easy to find solutions by hand:
diff(f(t),t,t)+omega^2*f(t)=0.
However, I solve numerically.

restart;
ode:=diff(f(t),t,t)+omega^2*f(t)=0;
res:=dsolve({ode,f(0)=f0,D(f)(0)=0},numeric,output=listprocedure,parameters=[f0,omega]);
F,F1:=op(subs(res,[f(t),diff(f(t),t)]));
p:=proc(a,b,f0,omega) local fa,f1a,fb,f1b,i; global p1,p2,p3,p4;
   res(parameters=[f0,omega]);
   p1(_passed):=F(a);
   p2(_passed):=F1(a)+1;
   p3(_passed):=F(b);
   p4(_passed):=F1(b)-1;
   if type(procname,indexed) then cat(op(0,procname),op(procname))(_passed) else [p1,p2,p3,p4](_passed) end if
end proc;
assign(seq(cat(p,i)=subs(ii=i,(proc(a,b,f0,omega) p[ii](_passed) end proc)),i=1..4));
eval(p3); #Just checking the assignment
p(Pi/2,-Pi/2,1,1); #Checking p
p1(Pi/2,-Pi/2,1,1); #Checking p1
#Now solving f(a)=0, f'(a)=-1, f(b)=0, f'(b)=1 for a,b,f0,omega:
fsolve([p1,p2,p3,p4],[1.5,-1.5,1,1]);
res(parameters);
plots:-odeplot(res,[[t,f(t)],[t,diff(f(t),t)]],-3..3); #One solution
fsolve([p1,p2,p3,p4],[3,-3,.5,1]);
plots:-odeplot(res,[[t,f(t)],[t,diff(f(t),t)]],-3..3); #Another solution

It seems that there is a diff versus Diff problem. The number option is only used when the system is given by a procedure.
Try this instead:

DEplot3d(value(sys), [u1(t), u2(t), u3(t)], t = 0 ..0.4,[[u1(0) = 1, u2(0) = 1, u3(0) = 3]]);
# and
DEplot(value(sys), [u1(t), u2(t), u3(t)], t = 0 ..0.4,[[u1(0) = 1, u2(0) = 1, u3(0) = 3]],scene=[u1(t),u2(t)]);

Will the following device serve your purpose?

u:=sqrt(2)*sqrt(g)/sqrt(b);
evalindets(u,posint,``);
combine(%) assuming positive; #One square root
expand(%); #Back to positive integers factored out

Addressing the problem you mention at the end, here is what I did.

Q:=proc(pp2,fi0) local res,F0,F1,F2,a;
global Q1,Q2;
print(pp2,fi0);
if not type([pp2,fi0],list(numeric)) then return 'procname(_passed)' end if:
res := dsolve(subs(p2=pp2,phi0=fi0,{eq1=0,eq2=0,eq3=0,u(0)=0,u(1-zet)=0,phi(0)=phi0,T(0)=0,D(T)(0)=1}), numeric,output=listprocedure):
F0,F1,F2:=op(subs(res,[u(eta),phi(eta),T(eta)])):
a[1]:=evalf(2/(1-zet^2)*Int((1-eta)*(F1(eta)*rhop+(1-F1(eta))*rhobf)*( F1(eta)*rhop*cp+(1-F1(eta))*rhobf*cbf )/(F1(eta)*rhop+(1-F1(eta))*rhobf)*F0(eta),eta=0..1-zet))-pp2*10000;
a[2]:=evalf(2/(1-zet^2)*(Int((1-eta)*F0(eta)*F1(eta),eta=0..1-zet))) / evalf(2/(1-zet^2)*(Int((1-eta)*F0(eta),eta=0..1-zet)))-0.02;
Q1(_passed):=a[1];
Q2(_passed):=a[2];
if type(procname,indexed) then a[op(procname)] else a[1],a[2] end if
end proc;

Q1:=proc(pp2,fi0) Q[1](_passed) end proc;
Q2:=proc(pp2,fi0) Q[2](_passed) end proc;
fsolve([Q1,Q2],[6.5,.007]);
#On my not very powerful computer it took fsolve roughly 45 minutes to return the result:
[6.3329718919924559, 0.79154816988391156e-2];


Addendum.
Instead of using fsolve you may try to use minimization of the sum of squares of Q1 and Q2.
This can be done as follows:
Optimization:-LSSolve([Q1,Q2],initialpoint=[6.3,0.008]);
#The result agrees very well with the fsolve result.
#Now I did use a better initial point. But if I start with the same as in fsolve I get the same result in just about 2 minutes, i.e. more than 20 times as fast as fsolve:
Optimization:-LSSolve([Q1,Q2],initialpoint=[6.5,0.007]);
A somewhat speedier version uses the fact that you really need only compute 2 integrals not 3, since one of the integrals can be written as a linear combination of the other 2:
Q:=proc(pp2,fi0) local res,F0,F1,F2,a,INT0,INT10,B;
global Q1,Q2;
print(pp2,fi0);
if not type([pp2,fi0],list(numeric)) then return 'procname(_passed)' end if:
res := dsolve(subs(p2=pp2,phi0=fi0,{eq1=0,eq2=0,eq3=0,u(0)=0,u(1-zet)=0,phi(0)=phi0,T(0)=0,D(T)(0)=1}), numeric,output=listprocedure):
F0,F1,F2:=op(subs(res,[u(eta),phi(eta),T(eta)])):
INT0:=evalf(Int((1-eta)*F0(eta),eta=0..1-zet));
INT10:=evalf(Int((1-eta)*F0(eta)*F1(eta),eta=0..1-zet));
B:=(-cbf*rhobf+cp*rhop)*INT10+ rhobf*cbf*INT0;
a[1]:=2/(1-zet^2)*B-10000*pp2;
a[2]:=INT10/INT0-0.02;
Q1(_passed):=a[1];
Q2(_passed):=a[2];
if type(procname,indexed) then a[op(procname)] else a[1],a[2] end if
end proc;