In the first session l1 is renamed to the escaped local l1~. Then a=sqrt(l1~) is saved.

In the second session what is read in is a as sqrt(l1~) (actually (l1~)^(1/2)).

You can access that l1~ in (at least) two ways:

L:=op(1,a); #You could use l1 instead of L
subs(L=2,a); #Works

##  alternatively (and this time defining l1 as the old globalized escaped local):

l1:=cat(l1,`~`);
subs(l1=4,a);
simplify(%);

I don't see any evidence that assumptions are stored, though.

## OK one more way without doing any of the other is redefining a:
a:=subsop(1=l1,a);

##################################
Since you may have several assumed variables, you may wish a programmatic approach.
Example:

restart;
assume(l1>0,l2>0,l3>0); #Any assumptions
expr := sqrt(l1)+sin(x*l2)+exp(-l3)*l1;
save expr, "E:/test3.m";
#############
restart;
read "E:/test3.m";
expr;
nms:=convert(indets(expr,name),list);
nmsS:=convert~(nms,string);
L:=StringTools:-Substitute~(nmsS,"~",""); #Removing "~"
L1:=parse~(L);
S:=nms=~L1;
Expr:=subs(S,expr);
eval(Expr,{l1=4,l2=5,l3=6});
simplify(%);

In the help page for LinearAlgebra:-ConditionNumber we find the statement:

For floating-point Matrices, the norm of the inverse is computed by using estimates of the norms of the inverses of the triangular factors of the Matrix obtained by LU or Cholesky factorization.

Thus the method used is different from the one used in the non float case. After a few comparisons this difference doesn't seem to have anything to do with roundoff, but is rather a somewhat different measure of the condition number. Whether this should be considered a bug or just an irrelevant difference I shall leave to others to decide. (After all you might as well have used another norm).

In the exact (non-float) case this is used:

use LinearAlgebra in
Norm(A, infinity)*Norm(MatrixInverse(A), infinity)
end use;

If you repeat that with B you get 6.

If you try

lprint(u);

you find

21.1896000000000*z(t)-.784800000000000*sin(theta(t))*(0.5000000000e-1*2^(1/2)*cos(`thet1(t)`)+0.5000000000e-1*2^(1/2)*sin(`thet1(t)`))+.784800000000000*cos(theta(t))*sin(phi(t))*(-0.5000000000e-1*2^(1/2)*cos(`thet1(t)`)+0.5000000000e-1*2^(1/2)*sin(`thet1(t)`))-.784800000000000*sin(theta(t))*((1/2)*2^(1/2)*(.10*cos(`thet1(t)`)*cos(`thet2(t)`)-.10*sin(`thet1(t)`)*sin(`thet2(t)`)+.2*cos(`thet1(t)`))+(1/2)*2^(1/2)*(.10*sin(`thet1(t)`)*cos(`thet2(t)`)+.10*cos(`thet1(t)`)*sin(`thet2(t)`)+.2*sin(`thet1(t)`)))+.784800000000000*cos(theta(t))*sin(phi(t))*(-(1/2)*2^(1/2)*(.10*cos(`thet1(t)`)*cos(`thet2(t)`)-.10*sin(`thet1(t)`)*sin(`thet2(t)`)+.2*cos(`thet1(t)`))+(1/2)*2^(1/2)*(.10*sin(`thet1(t)`)*cos(`thet2(t)`)+.10*cos(`thet1(t)`)*sin(`thet2(t)`)+.2*sin(`thet1(t)`)))

You will notice that thet1(t) doen't appear as such, but as `thet1(t)`. This means that you have made it into a name probably unintentionally. In the present context, however, this is not so bad an idea. Then you don't even need Physics:-diff, you can use ordinary diff:

diff(u,`thet1(t)`);

You have done the same with thet2(t).

A simple way of changing the original u to the intended version:
w:=evalindets(u,name,parse);
lprint(%);
diff(w,thet1(t)); #Now this gives an error
Physics:-diff(w,thet1(t)); #This gives the result

I used piecewise in this. `dsolve/bvp/convertsys` doesn't seem to like that: It issues an error about not being able to convert to an explicit first order system. See note at the bottom.
After that I tried an ivp problem. dsolve had no problem with that.

restart;
sys:={diff(S(t), t) = -eta*K(t)*S(t)/(w*N*(S(t)+K(t))), diff(K(t), t) = eta*K(t)*S(t)/(w*N*(S(t)+K(t)))-piecewise(psi[S](t)>-c,upsilon,0),
diff(psi[S](t), t) = eta*K(t)^2*(psi[S](t)-psi[Iota](t))/(w*N*(S(t)+K(t))^2), diff(psi[Iota](t), t) = eta*S(t)^2*(psi[S](t)-psi[Iota](t))/(w*N*(S(t)+K(t))^2)};

bcs:={S(0)=1,K(0)=1,psi[S](T)=1,psi[Iota](T)=1};

params:={N=1,w=1,eta=1,upsilon=1,c=0.5,T=1};
res:=dsolve(eval(sys union bcs,params),numeric); #Error mentioned above
##The system looks good to me:
eval(sys union bcs,params);
# So try an ivp:
ics:={S(0)=1,K(0)=1,psi[S](0)=1,psi[Iota](0)=5};
res:=dsolve(eval(sys union ics,params),numeric); #No problem
plots:-odeplot(res,[t,psi[S](t)],0..2);
plots:-odeplot(res,[t,S(t)+K(t)],0..2.1,thickness=3);
plots:-odeplot(res,[t,S(t)],0..2);

NOTE: The problem `dsolve/bvp/convertsys` had was apparently caused by the use of S in psi[S]. S is thus appearing in two roles, as the S in S(t) and as the S in psi[S]. Changing the latter to e.g. s (lower case s) makes dsolve work on the bvp.

I shall consider y1 just a third variable until I hear otherwise (see my comment above).

Furthermore I will not be using the name D for the operator, since it is used by Maple already.

Q:=f->unapply(convert(diff(f(x,y,y1),x)+y1*diff(f(x,y,y1),y)+y*y1/x*diff(f(x,y,y1),y1),D),x,y,y1);

Example:

L:=(x,y,y1)->x^2*sin(y)+37*y1;

Q(L)(s,t,u);

Q(L)(8,Pi/2,-6);

## If you prefer you can use a more explicit procedural form, which also has the advantage that x,y,y1 can be declared local to Q:

Q:=proc(f) local expr,x,y,y1;
  expr:=diff(f(x,y,y1),x)+y1*diff(f(x,y,y1),y)+y*y1/x*diff(f(x,y,y1),y1);
  expr:=convert(expr,D);
  unapply(expr,x,y,y1);
end proc;

### Comment:
One could have expected the more direct definition
Q:=f->unapply(D[1](f)(x,y,y1)+y1*D[2](f)(x,y,y1)+y*y1/x*D[3](f)(x,y,y1),x,y,y1);

would work, but it doesn't. Thus the extra step via diff.

You have to realize that (-1)^(1/3) is exp(I*Pi/3)=1/2+1/2*I*sqrt(3), i.e. it is the principal value of the cube root of -1. Thus it is imaginary.
To get the realvalued cube root use surd:

plot(surd(x,3),x=-1..1);
#Or just use convert:
plot(convert(x^(4/3)*sin(1/x),surd),x=-1..1);

See
?surd

This is taking from the example you gave in

http://www.mapleprimes.com/questions/202124-Defining-A-Function-Of-A-PreDefined#comment222568

restart;
n:=7:
g_temp := add((t[i]+1)*ln(t[i]+1)-t[i]*ln(t[i]), i = 1 .. n);
g := unapply(g_temp, seq(t[i], i = 1 .. n));
##First version is cleaner than the second (in my opinion):
p:=proc(x::Vector) local N; global g,step;
  N:=op(1,x);
  x+step*Vector(N,i->D[i](g)(seq(x[j],j=1..N)))
end proc;
p(<1,2,3,4,5,6,7>);
## This uses the arrow notation:
NS:=(x::Vector)->x+step*Vector(op(1,x),i->D[i](g)(seq(x[j],j=1..op(1,x))));
NS(<1,2,3,4,5,6,7>);

You could use this procedure mp instead of modp:

mp:=(x,p)->x-floor((x/p))*p;

odeplot(soln, [mp(theta(t), 2*Pi), omega(t)/`&omega;H`], 0 .. n*Th, labels = ["&theta;(t)","&omega;(t)/&omega;H"], axes = boxed, style = point, size = [.25, .75])

Another version using the builtin trunc:

mpt:=(x,p)->x-`if`(x>=0,trunc(x/p),trunc(x/p)-1)*p;

I vaguely remember a discussion about this on MaplePrimes concerning which version of an extension of modp to reals is better/faster.

Certainly you cannot have an inequality as an initial condition. The error message actually says that, but in a positive way: it tells you what you can have, i.e equations or expressions. Giving an expression 'expr', say, will be interpreted as giving the equation expr=0.

Does this do it:

display(seq(HyperionOrbit(i, `&omega;H`*i, 1, 1), i = -3 .. 3), view = [-3 .. 3, default]);

It is OK in Maple 18.
A look at the code in Maple 18 and Maple 2015 shows that the code has been changed a bit.
In Maple 2015 there is a puzzling line 4:

hor := evalb(`if`(orientation <> (':-default'),orientation = (':-horizontal'),horizontal = true and vertical = false));

which seems to me to return false if orientation=default, because vertical is per default true.

The value of hor is used at the end of the code

  88     if hor then
  89         vp[i] := [-y, x]
           else
  90         vp[i] := [x, y]
           end if;

A workaround is to give the orientation:

DrawNetwork(N,orientation=horizontal);

##Another one (not so nice).
DrawNetwork(N,horizontal=true,vertical=false);
DrawNetwork(N,horizontal,vertical=false); #Shorter, but slightly confusing.

The code would have worked as intended if the default for 'vertical' had been vertical::truefalse := not horizontal

I corrected your equations some and used numerical solution with the parameters option:

restart;
xm:=t*vm; #Position of master
ode1:=diff(f1(t),t)=(xm-f1(t))/sqrt((xm-f1(t))^2+f2(t)^2)*vd;
ode2:=diff(f2(t),t)=-f2(t)/sqrt((xm-f1(t))^2+f2(t)^2)*vd;
res:=dsolve({ode1,ode2,f1(0)=0,f2(0)=h},numeric,parameters=[h,vm,vd]);
res(parameters=[10,1,5]); #Setting the parameters
plots:-odeplot(res,[f1(t),f2(t)],0..2.1);
plots:-odeplot(res,[f1(t),f2(t)],0..2.1,view=[0..2.1,0..10],frames=25); #An animation: Click and run

Note: When the dog joins its master they have the same position, so the square root is zero. This creates the numerical problem at t = 2.0834...

###  The dog follows along after it joins its master (using piecewise):
restart;
xm:=t*vm;
ode1:=diff(f1(t),t)=piecewise(f2(t)>0,(xm-f1(t))/sqrt((xm-f1(t))^2+f2(t)^2)*vd,vm);
ode2:=diff(f2(t),t)=piecewise(f2(t)>0,-f2(t)/sqrt((xm-f1(t))^2+f2(t)^2)*vd,0);
res:=dsolve({ode1,ode2,f1(0)=0,f2(0)=h},numeric,parameters=[h,vm,vd]);
res(parameters=[10,1,5]);
plots:-odeplot(res,[f1(t),f2(t)],0..3,thickness=3);
plots:-odeplot(res,[f1(t),f2(t)],0..5,view=[0..5,0..10],thickness=3,frames=25);

### Although not a dog owner I tried a third version using events instead of piecewise. Also T(t) is given the time of the happy occasion/collision:
restart;
xm:=t*vm;
ode1:=diff(f1(t),t)=(xm-f1(t))/sqrt((xm-f1(t))^2+f2(t)^2)*vd*b(t)+(1-b(t))*vm;
ode2:=diff(f2(t),t)=-f2(t)/sqrt((xm-f1(t))^2+f2(t)^2)*vd*b(t);
res:=dsolve({ode1,ode2,f1(0)=0,f2(0)=h,b(0)=1,T(0)=0},numeric,discrete_variables=[b(t)::boolean,T(t)::float],parameters=[h,vm,vd],events=[[f2(t)=0,[b(t)=0,T(t)=t]]]);
res(parameters=[10,1,5]);
plots:-odeplot(res,[f1(t),f2(t)],0..3,thickness=3);
res(3);
plots:-odeplot(res,[f1(t),f2(t)],0..5,view=[0..5,0..10],thickness=3,frames=25);

Replace the assignment to u(k) with an assignment to some other local variable, e.g. v:

v:=max(exp((0.5*(k_r+1)^2)*0)*(exp(0.5*n*dx*(k_r+1))-exp(0.5*n*dx*(k_r-1))),0);
u:=[op(u),v];

Your u is not an array, it is a list. If you do want an array you could do:

am:=proc(dx,Nminus,Nplus)
local k_r,n,a2,u;
k_r:=1/2;
a2:=a/2;
u:=Array(Nminus..Nplus);
for n from Nminus to Nplus do
u[n]:=max(exp((0.5*(k_r+1)^2)*0)*(exp(0.5*n*dx*(k_r+1))-exp(0.5*n*dx*(k_r-1))),0);
end do;
u
end proc;

Or with indexing from 1:

am:=proc(dx,Nminus,Nplus)
local k_r,n,a2,u;
k_r:=1/2;
a2:=a/2;
u:=Array(1..Nplus-Nminus+1);
for n from Nminus to Nplus do
u[n-Nminus+1]:=max(exp((0.5*(k_r+1)^2)*0)*(exp(0.5*n*dx*(k_r+1))-exp(0.5*n*dx*(k_r-1))),0);
end do;
u
end proc;

You are not using numerical integration. You are computing the symbolic answer (using int) and then you are applying evalf.
Try replacing int by Int. Then you will have numerical integration.
Also try removing evalf, but using int.

Thirdly, try

int((exp(x)*(4420*cos(4)*sin(4)-544*cos(4)^2+148147*exp(-1)-4225*cos(4)-215203)/(71825*exp(1)-71825*exp(-1))-exp(-x)*(4420*cos(4)*sin(4)-544*cos(4)^2+148147*exp(1)-4225*cos(4)-215203)/(71825*exp(1)-71825*exp(-1))+(32/4225)*cos(4*x)^2+(1/71825)*(4225+(2210*x-6630)*sin(4*x))*cos(4*x)+x^2+8434/4225)^2, x = 0 .. 1,method=_RETURNVERBOSE); ## I learned that option from acer.

(The only successful methods FTOC and FTOCMS are more easily seen by wrapping the whole thing in evalf).

The order and/or the arrangement of the terms in the exact result will influence a subsequent use of evalf, i.e. the rounding errors will differ.

In Eqv[i] there is a missing (t) after v[i] in Av*(1-r[i]^2)*v[i].

The parameter omega[st] is given no value.

The initial conditions for the v's are all zero. It follows from the odes for the v's that they will remain zero (by uniqueness).

I tried:
ic:={u[1](0)=0.8, v[1](0)=0.1,u[2](0)=0.8, v[2](0)=0.1,u[3](0)=0.8, v[3](0)=0.1,u[4](0)=0.8, v[4](0)=0.1};
and used omega[st]=50.

sys:=map(op,eval({seq(EqSys[i],i=1..4)},{params}));
res:=dsolve(eval(sys,omega[st]=50) union ic,numeric);
plots:-odeplot(res,[u[1](t),v[1](t)],0..1);

plots:-odeplot(res,[seq([t,v[i](t)+i/10],i=1..4)],0..1,thickness=2);