I've noticed that many of your Questions involve using solve in a for loop to solve for a variable that is absolutely trivial to solve for by hand. The following produces exactly the same results as Tom Leslie's version of your code, but it doesn't use solve:
 

Download NoSolve_(1).mw

Your procedure can be made faster by

1. using diff instead of D,
2. simplifying binomial(n,k)/n! to 1/((n-k)!*k!),
3. remembering and re-using the derivatives and factorials of smaller orders.

restart:
Taylor2V:= proc(f::algebraic, V::[name= algebraic, name= algebraic], N::posint:= Order)
local 
    v:= lhs~(V), V0:= (lhs-rhs)~(V)[], n, k, 
    `!`:= proc(n) option remember; n*thisproc(n-1) end proc,
    D:= proc(n,k)
    option remember;
        `if`(k=0, diff(thisproc(n-1, 0), v[1]), diff(thisproc(n, k-1), v[2]))
    end proc
;
    D(0,0):= f; `!`(0):= 1; #Initialize remember tables.
    add(add(eval(D(n-k,k), V)*`*`(V0^~(n-k, k))/`!`(n-k)/`!`(k), k= 0..n), n= 0..N)
end proc
:
f:= 1/(2+x*y^2):
S:= CodeTools:-Usage(Taylor2V(f, [x=0, y=0], 30)):
memory used=5.49MiB, alloc change=32.00MiB, 
cpu time=31.00ms, real time=32.00ms, gc time=15.62ms

This Answer is intended for your general education into how to make Maple procedures faster. But you'll never make this as fast as mtaylor because mtaylor (which is written in Maple) calls series (which is NOT written in Maple) using VV's trick.

The ssystem command returns a list of two items. The first item is an error code (0 meaning no error); the second item is a string containing the program's output. This string is the only part that you want.

The commands Import and ImportGraph have different syntax, and you are mixing them up. Import will take raw string data; ImportGraph requires a file or filename. Try this:

L2:= Import(
    ssystem("D:/nauty27r3/geng -c -b 6 -g")[2], 
    source= direct, format= "Graph6", output= list
);

Note that the command does not begin with "!".

However, you can verify that C6 is embeddable in S5 with a single command:

GT:= GroupTheory:
GT:-MinPermRepDegree(GT:-CyclicGroup(6));
                               5
 

Here's a procedure for it:
 

Download PolyToList.mw

I hope this helps you:
 

Download OrbitPlot.mw

Like this:

z:= [
    0, 0, 0, 0, 0, 0, .689376362, 1.378752724, 2.068129087, 2.757505449, 0, 1.02920355,
    2.0584071, 3.087610649, 4.116814199, 0, 1.216469264, 2.432938529, 3.649407793,
    4.865877057, 0, 1.325720912, 2.651441823, 3.977162735, 5.302883646
]:
plots:-listcontplot(
    [ListTools:-LengthSplit](z,5), axis[1,2]= [tickmarks= ([$1..5]=~[$0..4])]
);

 

The Maple and Mathematica solutions that you show are mathematically equivalent. Both will give nonreal values for any real values of the a's and b's (unless a[1,1]*b[0,2] = a[0,2]*b[1,1], in which case the discriminant is 0). The expression under the square root sign cannot be positive for any real values of the parameters because it factors as (-1)*(a[1,1]*b[0,2] - a[0,2]*b[1,1])^2.

The absence of I from an expression is not a reliable indicator that the expression is real-valued. The presence of I is not a reliable indicator that it is not real-valued.

For rational-coefficient polynomial equations (such as yours) with parameters, you can get solve to break down all the relevant real intervals of the parameters which produce real solutions by adding the options parametric and real:

solve(
    k^2*a[1,1]^2 + k^2*b[1,1]^2 + 4*k*a[0,2]*a[1,1] + 
    4*k*b[0,2]*b[1,1] + 4*a[0,2]^2 + 4*b[0, 2]^2 = 0, 
    k, parametric, real
);

Note that the solution below covers all degenerate cases, i.e., when the equation is not really quadratic because the coefficient of k^2 is 0 and even when the coefficient of k is also 0. The parametric option always considers those cases.
   

I agree that a list or sequence of vectors would be the preferred format for programmatic usage. For a displayable table, I prefer using a matrix. (I especially like the way that matricies display with their entries centered in both dimensions.) The code below does both: My procedure DiffTab returns the "raw" list of vectors, the same as Tom Leslie's doDiff. The procedures PadV and DisplayTab turn those vectors into a neatly displayable matrix.

PadV:= proc(V::Vector, n::posint)
local k:= numelems(V), OP:= j-> n-k-1+2*j;
    Vector(2*n-1, applyop~(OP, 1, op(2,V)), fill= ``)
end proc
:
DisplayTab:= (DT::list(Vector), X::Vector)-> `<|>`(PadV~([X, DT[]], numelems(X))[])
:
DiffTab:= proc(Y::Vector)
local k, r:= Y, R:= table([1=r]);
    for k from 2 to numelems(Y) do r:= r[2..] - r[..-2];  R[k]:= r od;
    convert(R, list)
end proc
:
X:= <[$1..5]>:
Y:= 1/~X:
DT:= DiffTab(Y);
<<` X` | ` Y` | [``$4][]>, DisplayTab(DT,X)>;

You need to change xy to x*y.

Note the diffferent output from these two:

restart:
solve(y(1)^2 = 2., y(1));
             1.414213562, -1.414213562

{ fsolve(y(1)^2 = 2., y(1)) };
             { }

The fsolve returns no output (which I emphasized by putting the command in { }), not even an error message. Technically, the fsolve output is called NULL.

An expression like y(1) is called in Maple an unevaluated function call (or sometimes just a "function"). So, you can use a function as a variable for solve but not for fsolve. So, your call T__v_x(1000) is producing NULL output in a position where some other command is expecting a number.

Do you mean something like this?

#x and y coordinates of grid:
GX:= [seq](0..3, 1/3):  GY:= [seq](1..3, 1/3)
:
#x and y coordinates of highlighted points:
PX:= [seq](1/4..11/4, 1/4):  PY:= [seq](5/4..11/4, 1/4)
:
plots:-display(
    #border:
    plot(
        [[GX[1],GY[1]], [GX[1],GY[-1]], [GX[-1],GY[-1]], [GX[-1],GY[1]], [GX[1],GY[1]]],    
        color= COLOR(RGB, .3$3), thickness= 2
    ),
    #inner:
    plot(
        [
            seq([[GX[1],y],[GX[-1],y]], y= GY[2..-2]), 
            seq([[x,GY[1]],[x,GY[-1]]], x= GX[2..-2])
        ],
        color= COLOR(RGB, .5$3), thickness= 0
    ),
    #points:
    plot(
        [seq](seq([x,y], x= PX), y= PY),
        style= point, symbol= solidcircle, symbolsize= 8,
        color= red
    ),
    axes= none, scaling= constrained
);

Your general solution isn't general enough. If f is *any* differentiable function such that f(0)=1, then x0 = y0 = f is a solution to your equation. 

A basic fact about systems of equations is that the number of things that can be solved for is at most the number of equations. 

If L is the list of expressions, then all that you need to do is

R:= Grid:-Map(int, L, x= 0..1);

For numeric integration, do 

R:= Grid:-Map(int, L, x= 0..1, numeric);

The command will automatically apply some rudimentary load balancing.

It's possible to get the exact solutions for all 6 functions by direct use (i.e., without using dsolve) of the Method of Undetermined Coefficients. If you try to do this with dsolve, it will get the first 5, but the 6th took longer than I cared to wait. I killed it after more than an hour. The method below gets it in seconds.
 

Download UndeterminedCoeffs.mw