VV's Answer shows you how to get the generalized eigenvectors that you want.

The purpose of this Answer is to correct several of your misconceptions about defective matrices.

You wrote:

• When a real matrix have repeated eigenvalue (i.e. multiplicity >1) and the matrix is not the identity matrix, then it is called defective eigenvalue. 

A defective eigenvalue is one for which the space spanned by the corresponding eigenvectors has dimension less than the (algebraic) multiplicity of the eigenvalue. While it's true the such eigenvalues always have multiplicity greater than 1, the converse is not true: A repeated eigenvalue is not necessarily defective. There's nothing special about the identity matrix in this regard. There's nothing special about real matrices.

• There is an algorithm to generate another (linearly independent) eigenvector from this same eigenvalue called the defective eigenvalue method.

The algorithm that you speak of does not generate (true) eigenvectors; rather, it generates generalized eigenvectors. Given the matrix A and vector v__2 from your example, you can easily check that A.v__2 <> 3*v__2, and thus v__2 is not an eigenvector.

Maple's Eigenvectors does return a maximal linearly independent set of true eigenvectors. 

I assume that the data in the Excel file is a two-column matrix, the first column being times and the second being values of Q. If that's not the format, slight adjustments to the following will be needed.

First, import the Excel file to a Maple matrix:

QX:= ExcelTools:-Import(filename):

Then,

QI:= Interpolation:-LinearInterpolation(QX);
Q:= t-> `if`(t::numeric, QI(t), 'procname'(t)):

Set up your ode system odesys just like you have above, including the Q(t). Then

sol:= dsolve(odesys, numeric, known= Q);

The process can be easily generalized to exchange any "outer" function or set of functions with any "inner" function or set of functions, and to do it at arbitrary functional positions:

ExchangeFuncs:= proc(
    e, 
    F::[{name, set(name)}, {name, set(name)}],
    p::[posint, posint]:= [1,1]
)
    evalindets(
        e,
        And(
            specfunc(anything, F[1]), 
            'patfunc'(
                anything$(p[1]-1), 
                And(specfunc(anything, F[2]), 'patfunc'(anything$p[2])), 
                anything
            )
        ),
        u-> 
            op([p[1],0], u)(
                op([p[1], ..p[2]-1], u),
                op(0,u)(op~([..p[1]-1, p, p[1]+1..], u)[]), 
                op([p[1], p[2]+1..], u)
            )
    )
end proc
:    
#Example (inner functions are outer functions' 2nd argument):
e:= f(1, g(2,3)) + h(4, i(5,6)):
#Apply outer functions to inner functions' 1st argument:
ExchangeFuncs(e, [{f,h}, {g,i}], [2,1]); 
                 g(f(1, 2), 3) + i(h(4, 5), 6)

#Apply outer functions to inner functions' 2nd argument:
ExchangeFuncs(e, [{f,h}, {g,i}], [2,2]);
                 g(2, f(1, 3)) + i(5, h(4, 6))

To do the job that you requested:

MySumOfInt:= ExchangeFuncs(MyIntOfSum, [Int, Sum]);

Or you may want to make the second argument [{Int, int}, {Sum, sum}].

You can reverse the exchange, which'll return the original:

ExchangeFuncs(MySumOfInt, [{Sum, sum}, {Int, int}]);

There's no need to do something different in the code itself. The plotting can be done as an afterthought:

plots:-listplot([seq(A[i], i= 1..1000)]);

Use a digraph in GraphTheory, like this:

GT:= GraphTheory:
G:= GT:-Graph({
    [[1,2], 35], [[3,4], 20], [[6,5], 15], [[8,7], 40],
    [[2,8], 222], [[4,2], 111], [[6,4], 333], [[8,6], 444]
}):
GT:-SetVertexPositions(
    G,
    [[-1,2], [0,1], [2,2], [1,1], [2,-1], [1,0], [-1,-1], [0,0]]
):
plots:-display(
    subs(
        ["111"= x__1, "222"= x__2, "333"= x__3, "444"= x__4], 
        GT:-DrawGraph(
            G, 
            stylesheet= [
                edgefont= [times,bold,16],
                vertexcolor= black,
                vertexfontcolor= black,
                vertexfontsize= 6
            ]
        )
    ),
    view= [(-.75..1.75)$2]
);

All sizes, colors, and positions can be adjusted, of course.

Okay, so you mean that you want a vector-field plot of the gradient. That's do-able, and shown below. To understand my comment about streams, you need to look at that procedure by expanding the Code Edit Region in the worksheet that you posted. Here it is:

streams:=proc(n,fcn,pts,name) 
       local f,k,DISK,p; 
       f:=transform((r,th) -> [r*cos(th),r*sin(th)]): 
       DISK:=disk([0,0],1,color=red); 
       for k from -n by 2 to n do 
             p[k]:=implicitplot(fcn=2*k/n,x=-n/2..n/2,y=-Pi..Pi,numpoints=pts): 
       od: 
       display(seq(f(p[2*i]), i=-n/2..n/2), DISK, view=[-n/2..n/2,-n/2..n/2], scaling=constrained, axes=boxed, title=name); 
end:

As you can see, it uses implicitplot, which doesn't work for vector fields.

To do the plot of the gradient, first take out (or comment out) these lines:

with(Physics[Vectors]):
Setup(mathematicalnotation= true):
V:= Gradient(V):

And add these lines:

dV:= VectorCalculus:-Gradient(V, [x,y]):
plots:-fieldplot(dV, x= -2..2, y= -2..2, fieldstrength= log[10]);

Then re-execute the worksheet (so that the restart is used). You'll get this field plot:

There are many options to control the appearance, size, and relative size of the arrows. See the help page ?plots,fieldplot for details.

I don't know what you meant by scene= [P, P(t)]. I changed it to scene= [t, P(t)]. I don't know what you meant by conditions. I just omitted that. If the below is not what you want, then if you can explain what you meant, I'll probably be able to tell you how to do it.

DEtools:-DEplot(
   diff(P(t),t) = P(t)*(7-P(t))-10, P(t), t= -10..10, P= -10..10, 
   scene= [t, P(t)], number= 1, stepsize= .01
);

To include a trajectory in the plot, you of course need intial conditions,

Let's say that LL is a list of lists of real numbers; it doesn't matter how many. Then do

plot(map(ListTools:-Enumerate, LL));

Presumably, ff is a procedure, and we'd need to see its full code to determine its "thread safety". If you've experienced weird behavior with this, a lack of thread safety is the likely cause. 

You can implement that formula with overload, like this:

restart:
`Psi/original`:= eval(Psi):
unprotect(Psi):
Psi:= overload([
    proc(r::And(positive, fraction), $)
    option overload;
    local n:= trunc(r), f:= r-n, p, q, k;
        (p,q):= (numer,denom)(f);
        simplify(
            q*add(`/`~([seq](p..(n-1)*(q+1), q))) +
            2*add(cos(2*Pi*k*f)*ln(sin(Pi*k/q)), k= 1..iquo(q-1,2)) - 
            Pi/2*cot(Pi*f) - ln(2) - ln(q) - :-gamma
        )
    end proc,

    `Psi/original`
]):
protect(Psi, `Psi/original`)
:
Psi(1/12);

Or, if you want more control over when the conversions are performed, you could do the exact same thing to `expand/Psi` rather than Psi.

Apply the expand command to the factored polynomial, the denominator in this case.

I think that you need to have some of Acer's points stated more emphatically: 

1. Your command print(pout) is irrelevant. All that it does is print the filename. 

2. Using the character terminating commands and loops (semicolon or colon) to control the plot output is error prone. You shouldn't use this method.

3. The command that actually produces the final plot that you want saved, display in this case, should be wrapped with a print command. 

1. There's an error in your input of the ODEs: The Pr(...) should be Pr*(...). Maple doesn't catch this type of error when Pr is assigned a numeric value (as it has been); unfortunately, this error simply causes the part in parentheses to be ignored. You have the same problem with alpha(...).

This code is obviously a variation of code originally posted by me. The numerous superfluous parentheses around derivatives indicate that it had been converted to 2D Input and then back to plaintext. I wonder how many hands this code has passed through containing these errors? 

2. You can vary the linestyle using syntax analogous to that used to vary the color. The valid linestyle names are dash, dashdot, dot, longdash, solid, spacedash, spacedot, and ticks. So, for example, linestyle= [solid, dash, dot].

3. This boundary-layer problem uses a finite value as a substitute for an upper boundary at infinity. This should be indicated by a caption on the plots. I think that your 1 is a very bad choice for this value because it obfuscates that the value is a fake upper boundary.

You should not specify a range for eta in the plotting commands. Since this is a BVP, the range will be automatically inherited from the dsolve command.

You should vary the value substituted for infinity to make sure that this does not cause too large a difference in the solutions. The validity of the results is suspect until this has been done. I've read many journal papers on boundary-layer problems that either ignore this point or fail to mention it. Note that the best choice for this value depends on the parameters. This is almost always ignored in these papers.

Maple's mesh-based BVP solvers often become unstable when too large a value is used for a fake infinity in a boundary-layer problem. The solvers can usually detect this instability, but they can't compensate for it; instead they usually return an error message. This is the most difficult issue with solving these boundary-layer problems with Maple.

I'm guessing that the thing that you find most confusing about the procedure is that it doesn't declare its parameters. Instead, it refers to them as args[1], args[2], and args[3]. The args and nargs are keywords: The former is the argument sequence, and the latter is the count of arguments.

Here is my version of the procedure with modern parameter specifications:

BasisGrid:= proc(
    u::{Vector(2, numeric), [numeric, numeric]},
    v::{Vector(2, numeric), [numeric, numeric]},
    ab::{Vector(2, numeric), [numeric, numeric]},
    {pointoptions::list(name, name= anything):=
        [symbol= circle, symbolsize= 20, color= blue]
    }
)
    local combo:= [seq](ab[1]*u+ab[2]*v);
    plots:-display(
        Lamp:-Gridvects([seq](u), [seq](v)),
        Lamp:-Gridvects(
            [1,0], [0,1], max(abs~(combo), 5),
            'vectors'= 'off'
        ),
        plot(combo, 'style'= 'point', pointoptions[]),
        _rest
    )
end proc
:

The _rest is a keyword for arguments that were passed but not declared.

for a to 2 do k||a:= a od;