I cannot fault Tom Leslie for having given up trying to understand what you were trying to do; your number of errors is extreme. However, having read through your code numerous times over the past 17 hours, I figured out that you are trying to find a degree-12 power series solution to the IVP

   y'''(x) = -1/2*y(x)*y''(x), y(0) = 0, y'(0) = 1, y''(0) = A.

That can be done like this:

Order:= 13:  # "environment" variable for degree of "order" term of power series 
S:= dsolve(
   {diff(y(x),x$3) = -1/2*y(x)*diff(y(x),x$2), y(0)=0, D(y)(0)=1, (D@@2)(y)(0)=A},
   y(x), 'series'
);
G:= unapply(convert(S, 'polynom'), x); 

Pay close attention to the where I have placed the parentheses! Also, parentheses and square brackets are not interchangable. On the other hand, the spacing and line breaks are a matter of personal style and can be changed however you desire.

Depending on the size of B and the uniqueness of its elements, there can be a considerable benefit to converting B to a set first:

remove(member, A, {seq(B)})

Let n be the number of elements of A, let m be the number of elements of B, and let q be the number of distinct elements of B. Then (from theoretical considerations alone) the number of operations to perform remove(member, A, B) is proportional to n*m, and those to perform remove(member, A, {seq(B)}) is proportional to n*log(q)+m*log(m).

The lowercase e can only be used in explicitly specified floating-point constants. If you need to use a variable, then use 10. as a base, as shown by Acer.

The L^p norm of a function f (x) over an interval a..b is an integral---int(abs(f(x))^p, x= a..b)^(1/p)---for 1 <= p < infinity. The limit of this as p--> infinity is the "infinity norm", which (for continuous f) is also the maximum of the absolute value. So, Tom Leslie's Reply is incorrect about the L^2 norm that you asked about. He's just inadvertently recomputed the infnity norm (which can be clearly seen from his numeric ranswer).

Here's a very short procedure to compute it for any p, including infinity:

LpNorm:= (f::algebraic, xab::(name= range(complexcons)), p::positive)->
   `if`(p=infinity, numapprox:-infnorm(f, xab), int(abs(f)^p, xab, numeric)^(1/p))
:

Using it for your example:

Digits:= 15: #Good enough for most purposes.
exact:= x^4-x^3:
app:= -2.220446049250313e-16+.19676728660169784*x+4.6491983150099205*x^2
         -3.9010072748158082*x^3+2.3619056303875987*x^4-.33198038154796344*x^5:
map[3](LpNorm, exact-app, x= 0..2, [1, 2, infinity]); 
      [6.36243961193564, 5.55700549477000, 6.94938751138337]
 

You need to also save the whole solution (i.e., the direct output of your dsolve(..., numeric) command). If you save a snippet of code that refers to some other variable that you've defined, then you need to save that other variable also. Since X_br refers to solution, that is the case here. Note that you can save any number of variables together in a single file using a single save command. Here's how to do it:

restart
:
solution:= dsolve({diff(xbr(t),t$2) = -xbr(t), xbr(0)= 0, D(xbr)(0)= 1}, numeric):
X_br:= tau-> eval(xbr(':-t'), solution(tau)):
save X_br, solution, "/Users/carlj/desktop/NumSol.mpl"
:
#Now move to your plotting worksheet.
restart
:
read "/Users/carlj/desktop/NumSol.mpl":
plot(X_br, 0..7);

There are three important, significant differences between my code and yours:

1.  I do not use output= listprocedure ^[*1].
2. My X_br is a procedure; yours is an expression.
3. I include solution in the save command.

And there's one minor significant difference: I used the procedure-form plot command

plot(X_br, 0..7);

rather than the expression-form

plot(X_br(t), t= 0..7);

If you have some good reason for needing the expression form, that's still possible. Let me know.

The only difference that using a ".m" internal-format file makes is that Maple is a tiny bit more efficient working with those files. Sorry about my incorrect previous Answer regarding ".m" files.

^[*1]Actually, I never use output= listprocedure. I think that it's a shoddy programming practice, and there's actually nothing that's made easier by it.

You are right to be suspicious: The way that you propose to iteratively modify sets is quite inefficient. What you want to do can be done efficiently^[*1] by less than a quarter of a line of code:

{Ug3[], seq(seq(u..6, -6), u= Ug3)};

Note that this doesn't require the unwanted numbers to be removed because they aren't generated in the first place.

You asked if sets are efficient. Sets and lists are very efficient if you don't try to repeatedly modify them. If you need to repeatedly or iteratively modify the structure, use a Vector or table, and if you ultimately need a set or list, convert it to such when you're done with the modifications.

^[*1]Tom Leslie's one-liner is slightly more efficient than my code above.

Make the filename "filename.m" instead of "filename.mpl". The .m form is called "internal format". In addition to the code of the procedure, this format also puts in the file all the global variables needed to  execute the procedure correctly. 

[Edit: I was wrong about this. Ignore this Answer.]

It can be easily done using undocumented features of the Typesetting package. For example:

plot(x^2, x= -2..2, title= Typesetting:-mtext("My title", color= red));

Although these features are undocumented, they are stable and widely known because they're closely based on good ol' HTML.

The return value of print (and also of lprint, printf, and userinfo) is always NULL. That it displays information on your screen is only a side effect  (using the technical definition of that term) . If you map a function whose return value is NULL over a list,  of course the result is an empty list. So, the [ ] is the return value of print~ applied to any list.

If you don't want that [ ], you can use print~([a,b])[ ]. This is more versatile than VV's suggestion of ending with a colon because it's an expression rather than a statement; hence, it can embedded inside other expressions.

The equations are linear, so LinearAlgebra:-LinearSolve can be used. After you've entered the numeric data, the rest is a one-liner. The following works regardless of the number of equations:

restart:
N:= < 11, 23, 42, 82 >:
R:= < 16, 397, 67, 2809 >:
#The rest is automatic regardless of the number of equations: 
n:= numelems(N):
S:= [seq(a[n-j], j= 1..n)]=~ seq(LinearAlgebra:-LinearSolve(`<|>`(seq(N^~(n-j), j= 1..n)), R));
 
S := [a[3] = 1056677/24673210, a[2] = -238802589/49346420, 
      a[1] = 7777678199/49346420, a[0] = -29341339187/24673210]

I've indexed the a variables so that their index corresponds to the exponent in the polynomials. This is arbitrary; you can change it, or change a to any other name, without needing to change anything else.

I dislike using assign in this situation because having used it usually makes some future operation more difficult. But, if you wish, you can do assign(S) at this point.

Maple has the commands sprintf and, better yet, nprintf (just like sprintf, but creates variable names instead of strings) but, like assign, I wouldn't use them for this situation.

It seems like you don't expend much effort trying to find errors on your own, because you have so many, and some of them are so obvious (such as unbalanced quotation marks, or exchanging P and p (refering to your previous Question)). You just complain here that your code "don't work". Well, it's not like we here at MaplePrimes have some automatic way to find the errors. We just apply the logical rules by hand, as you could learn to do.

I think that this accomplishes what you want:

plots:-display(
   seq(
      seq(
         [
            plottools:-polygon(
               [[i,j],[i,j+1],[i+1,j+1],[i+1,j]], 
               color= `if`(j::odd, magenta, white)
            ), 
            plots:-textplot([i+.5, j+.5, sprintf("%d",i*j)])
         ][],
         i= 1..10
      ), 
      j= 1..10
   ),
   axes= none
);

If instead you want the numbers to increase as you go down, then change the j range from 1..10 to -1..-10, -1 and change i*j to abs(i*j) (or -i*j).

Your code executes without error for me. The file is plaintext, and is easily readable in most browsers or editors; so I can assure you that it contains exactly one matrix without any syntax issues. The first three positive roots of the determinant are

[0.269373966731148794754757341200462674642148338669373608296012774382885541088668300033662376394882171035784896969388875153559058382051590463209287590804e-1, 0.990674477020063147446983680621401355331598629110322989492668885141230642823201150297702439572643428277168426737360736443335299618942000792413816520383e-1, .101120616665296170767949963559947228275460392424861673345779263832172742662163268222684076145562969014514409071557800322470075818964100491423123404014]

These results agree with VV's above. His technique is more robust in terms of round-off error.

Obviously the number of cycles of length k is closely related to the number of walks of length k that end at their own starting point. Clearly, this number of walks is the trace of the k^th power of the adjacency matrix. It only remains to remove the walks that aren't cycles and to divide out the 2*k permutations of each cycle. This formula does that:

#Number of k-cycles for a given adjacency matrix A
kCycles:= (A::Matrix, n::posint, k::And(posint, Not({1,2})))->
   add(
      (-1)^(k-i)*binomial(n-i,n-k) * 
         add(LinearAlgebra:-Trace(A[S,S]^k), S= combinat:-choose(n,i)), 
      i= 2..k
   )/2/k
:

Continuing with the example of the octahedron graph:

A:= GraphTheory:-AdjacencyMatrix(GraphTheory:-SpecialGraphs:-OctahedronGraph()):
n:= 6: #number of vertices 
seq(kCycles(A, n, k), k= 3..6);
                         8, 15, 24, 16
add([%]); #total number of cycles
                               63

These results agree with VV's explicit list of the cycles of this graph.

My little procedure above will gain a huge efficiency improvement by using remember tables for the traces, but I don't have time right now. If someone else chooses to undertake this, recall that a remember table cannot be indexed correctly by a mutable structure like a Matrix. So, use S and k as the table indices.

Here's a procedure a little more general (because it works for an arbitrary number of vectors) than what you asked for:

InSpan:= proc(V::{{set,list}(Vector[column]), Vector[column], Matrix}, w::Vector[column])
   try LinearAlgebra:-LinearSolve(`if`(V::Matrix, V, `<|>`(`if`(V::Vector, V, V[]))), w)
   catch "inconsistent": return false
   end try;
   true
end proc:

So, the first argument can be a column vector, a matrix, or a set or list of column vectors; and the second argument is a column vector. Example:

InSpan({ <1,2>, <3,4> }, <5,6>);

Regarding the terminology: The span of a set of vectors is the vector space of all linear combinations of those vectors. 

Several points to mention:

1. The command with PDE(tools) is nonsense to Maple. You mean with(PDEtools);. However, I encourage to not use this unless we determine that it's necessary. (In particular, PDEtools can be used to obscure the functional dependencies, which might be nice for the display of these equations, but makes it impossible for me to understand them and help you further.)
2. Using unprotect on a name doesn't make it safe to use that name! Change unprotect(gamma); to local gamma;.
3. Square brackets must not be used in place of parentheses.
4. I need to see the functional dependencies is order to help you further. For example, if v__1 is a function of x and t then you must write v__1(x,t) wherever they are used (or you can use v__1(t,x), as long as you use the same variable order consistently throughout the problem.