#
# Create a temporary directory where I can store the
# Maple repository. I'm using the 'temp' directory
# on a windows machine.
#
# I'm using *full* paths to the directory, just to
# avoid any possible confusion/ambiguity.
#
# OP may want to change the definition of 'reposit'
# to something suitable for his/her machine
#
# If the directory already exists, remove it to ensure
# a nice "clean" start
#
  reposit := "C:/temp/reposit":
  if   FileTools:-Exists(reposit)
  then FileTools:-RemoveDirectory(reposit, recurse=true)
  fi;
#
# Now create the empty directory
#
  mkdir( reposit ):

#
# Create an empty repository called 'foo.mla' in the
# above directory.
#
# Again I'm using full paths, just to avoid any possible
# confusion/ambiguity
#
  march( 'create', cat( reposit, "/foo.mla"), 100 ):

#
# Check the contents of the repository "foo.mla". At this
# stage it should be empty!
#
  march('list', cat( reposit, "/foo.mla"));

#
# Define the OP's procedure and save it to the repository
# created above
#
  canSpell := proc(w)
                   option encrypted;
                   local blocks, i, j, word, found, N;
                   blocks := Array([{"B", "O"}, {"X", "K"}, {"D", "Q"}, {"C", "P"}, {"N", "A"},
                                    {"J", "W"}, {"H", "U"}, {"V", "I"}, {"A", "N"}, {"O", "B"},
                                    {"G", "T"}, {"R", "E"}, {"T", "G"}, {"Q", "D"}, {"F", "S"},
                                    {"E", "R"}, {"F", "S"}, {"L", "Y"}, {"P", "C"}, {"Z", "M"}]);
                   word := StringTools:-UpperCase(convert(w, string));
                   N := numelems(blocks);
                   for i to length(word) do
                       found := false;
                       for j to N do
                           if word[i] in blocks[j] then
                              blocks[j] := blocks[N];
                              N := N-1;
                              found := true;
                              break;
                           end if;
                       end do;
                       if not found then
                          return false;
                       end if;
                   end do;
                   return true;
             end proc:

  savelib(canSpell, cat( reposit, "/foo.mla")):

#
# Recheck the contents of the repository
#
# Everything *seems* to be working
#
  march('list', cat( reposit, "/foo.mla"));
  

  restart;

#
# After creating the repository, I attached it to an email
# and sent it to myself. I then saved the attachment to my
# "Desktop", defined below.
#
# The value of 'myDesk' below will need to be changed by the
# recipient of OP's email to correspond to where (s)he saved
# the email attachment
#
  myDesk:="C:/Users/TomLeslie/Desktop":

#
# Check my default library path
#
  libname;
#
# Temporarily redefine the libname to include the "foo.mla"
# library on my desktop
#
  libname:=libname, cat( myDesk, "/foo.mla");

#
# Check the contents of the added repository
#
  march('list', cat( myDesk, "/foo.mla") );

#
# Check encryption
#
  eval(canSpell);

#
# Run the OP's procedure
#
  canSpell("good");

#
# Define a couple of simple functions which
#
# 1) Convert from "phasor" representation to
#    complex numbers
# 2) Convert from complex representation to
#    "phasor" representation
#
  fromP:= (r, theta)-> r*exp(I*Pi*theta/180);
  toP:= x-> evalf([abs(x), (180/Pi)*argument(x)]);
#
# Define "phasors" (and convert to complex
# reprentation)
#
  U__LI:= fromP(230, 0):
  I__LI:= fromP(20, -30):
  Z__L:=  fromP(0.097, 7.2):
  I__N:=  fromP(40,-120):
#
# Perform calculation, and convert to "phasor"
# representation
#
   ans:= toP(U__LI - I__LI * Z__L + I__N * Z__L);

Error, (in pdsolve/numeric/plot) unable to compute solution for t>HFloat(0.0):
matrix is singular

  restart:

  f:=7*exp(-3*t)-(3*x):
  s:=0.05:
  T:=0.1:
  n:=round(T/s):
  t0:=0:
  x[0]:=6:
  for i from 0 by 1 to n-1 do
    #
    # Note the square brackets enclosing
    # the multiple substitutions
    #
    # the simplify() 'wrapper' is just so that
    # exp(0) gets evaluated to a floating point
    # which isn't really necessary at this stage
    #
      x[i+1]:=simplify(x[i]+s*subs( [x=x[0], t=t0+i*s], f)):
  od;

  restart;
  interface(version);

  with(ThermophysicalData);
  with(ThermophysicalData[CoolProp]);
  fluid := "R717"; T__in := 310; p__in := 11*10^5; v__in := 10;  A := 0.005;
  h__in := Property(massspecificenthalpy, T = T__in, P = p__in, fluid);
  rho__in := Property(density, T = T__in, P = p__in, fluid);
  m := A*v__in*rho__in;
  p__out := 2*10^5;
  eq1 := h__out = Property("massspecificenthalpy", "temperature" = T__out, "P" = p__out, fluid);
  eq2 := rho__out = Property("D", "temperature" = T__out, "P" = p__out, fluid);
  eq3 := h__in + v__in^2/2 = h__out + v__out^2/2;
  eq4 := m = A*v__out*rho__out;
  res := fsolve({eq1, eq2, eq3, eq4});
#
# Try adding plausible(?) constraints on the ranges
# of the variables
#
  res := fsolve( {eq1, eq2, eq3, eq4},
                 {T__out=200..500, h__out=0..5e06, v__out=0..100, rho__out=0..5}
               );

  restart;

#
# Define function
#
  doDiff:= (p,q,r)-> diff(p(q), q$r);
#
# A couple of examples
#
  doDiff( f, x, 3);
  doDiff( g, y, 2);
  doDiff( h, z, n);

#
# But is this any better that using
# D-operator notation?
#
  D[1$3](f)(x);
  D[1$2](g)(y);
  D[1$n](h)(z);
  convert(%%%, diff);
  convert(%%%, diff);
  convert(%%%, diff);

  restart;
  with(plots):
  delta:=1:
  inf:=500*delta+1;
  sol:= dsolve( [ 2*diff(y(x),x,x,x)+y(x)*diff(y(x),x,x)=0,
                  y(0)=0,
                  D(y)(inf)=1,
                  D(y)(0)=delta*D[1,1](y)(0)
                ],
                numeric,
                approxsoln=[y(x)=x],
                maxmesh=4096,
                method=bvp[midrich]
              ):
  odeplot( sol, [x, y(x)], x=0..1);
  odeplot( sol, [x, y(x)], x=0..inf)

  restart:
#
# Make sure enough digits are used for calculation
# but restrict how many are displayed, just for
# readibility
#
  Digits:=50:
  interface(displayprecision=4):
#
# Since system is 'linear' probably(?) better to use
# LinearAlgebra:-LinearSolve() rather than the more
# "straightforward" fsolve(). so load the relevant
# package
#
  with(LinearAlgebra):
#
# OP's system
#
  C[0]:= 3.19153824321146142351956847947*tau[1]-19.1492294592687685411174108768*tau[2]+111.703838512401149823184896781*tau[3]+3.19153824321146142351956847947*tau[4]-44.6815354049604599292739587124*tau[5]+622.349957426234977586315853494*tau[6]:
  C[1]:= 51.0646118913833827763130956714*tau[2]-612.775342696600593315757148056*tau[3]+51.0646118913833827763130956714*tau[5]-1429.80913295873471773676667880*tau[6]:
  C[2]:= -1.06073680388443795908856507616+3.19153824321146142351956847947*tau[1]+53.1609155734306093706448370717*tau[2]+1672.89412862088744108725223170*tau[3]+3.19153824321146142351956847947*tau[4]+27.6286096277389179824882892361*tau[5]+1026.57792701153122226218722129*tau[6]:
  C[3]:= -1.08847004231036963538035920033+3.19153824321146142351956847947*tau[1]+62.6399144226357196540662623767*tau[2]+2040.52109049201342887896297462*tau[3]+3.19153824321146142351956847947*tau[4]+37.1076084769440282659097145411*tau[5]+1242.54090729537544551915515930*tau[6]:
  C[4]:= -1.05523181556926815105314303389+3.19153824321146142351956847947*tau[1]+72.7671212023804312453829273862*tau[2]+2472.93216226733267613216245895*tau[3]+3.19153824321146142351956847947*tau[4]+47.2348152566887398572263795506*tau[5]+1512.91667059477930731128800348*tau[6]:
  C[5]:= -.922876006485286011069063957991+3.19153824321146142351956847947*tau[1]+82.9822841707707093164204255644*tau[2]+2971.36790137532483139495115633*tau[3]+3.19153824321146142351956847947*tau[4]+57.4499782250790179282638777288*tau[5]+1847.90980220852701343747673000*tau[6]:
#
# Convert to a matrix system and attempt a solution
#
  A,b:= GenerateMatrix
        ( [seq( C[j], j=0..5)],
          [seq( tau[j], j=1..6)]
        );
  LinearSolve(A,b);

Error, (in LinearAlgebra:-BackwardSubstitute) inconsistent system

#
# Idle curiosity - check the reduced row echelon form
# of the augmented matrix. Three observations
#
# 1. the fact that row 6 is all zero implies that the
#    system is underdetermined - so any solution would
#    require a free parameter, however
# 2. row 5 is equivalent to the statement 0=1, so no
#    solution is possible
# 3. The fact that all of the coefficients in the RRF
#    *seem* to be *exact* makes one wonder exactly how
#    the original equation system was generated!
#
  M:=ReducedRowEchelonForm(<A|b>);
#
# Regenerate the system of equations from the reduced
# row echelon form. Check the fifth and sixth entry
# in the output list, which confirms the statements
# (1) and (2) above
#
  GenerateEquations(M, [seq( tau[j], j=1..6)]);

  restart;
#
# Implement the sieve of Eratosthenes. This will
# return all prime numbers up to the supplied
# argument (reasonably efficiently)
#
# EG on my machine all primes up to 1000000 can
# be generated in less than 1sec
#
  sieve := proc( Nval::integer)
                 local vec, k;
                 description "The Sieve of Eratosthenes";
                 uses ArrayTools, LinearAlgebra:
               #
               # Initialise a vector with all ones
               #
                 vec:= Vector[row]( Nval, 1 ):
               #
               # set the first entry in the vector to zero
               #        
                 vec[1]:= 0:
               #
               # Set all even entries with index >=4, to zero
               # (eliminates all the even indices which can't
               # be prime)
               #
                 Fill(0, vec, 3, 2);
               #
               # For each odd number, say k,  eliminate all
               # of its odd multiples (other than 1*k)
               #
                 seq( Fill
                      ( 0, vec, k^2-1, k ),
                      k=3..floor( sqrt(Nval) ),
                      2
                    ):
               #
               # return the indices for all non-zero values
               #
                 return Transpose
                        ( SearchArray
                          ( vec )
                        );
           end proc:
#
# Generate all the primes between 1 and N, Change
# the following assignment as necessary
#
  N:=1000000:
  prims:=sieve(N):
#
# Generate a random number between 1 and the
# number of primes in the vector 'prims'
#
  randomize():
  r:= rand( 1..op(1, prims)):

#
# Now generate (say) a hundred "random" prime
# numbers between 1 and N (defined above)
#
  seq( prims[ r() ], j=1..100);

delta:=-5:
inf:=1000;
ode:=2*diff(y(x),x,x,x)+y(x)*diff(y(x),x,x)=0:
ics:= y(0)=0, D(y)(inf)=1, D(y)(0)=delta*D[1,1](y)(0):
sol:=dsolve( [ode, ics], numeric, approxsoln=[y(x)=x], maxmesh=4096,method=bvp[midrich]);
odeplot( sol, [x, y(x)], x=0..1);
odeplot( sol, [x, y(x)], x=0..inf)