What's happening is that the copy-and-paste operation is overwhelming Maple's GUI component because it's trying to construct a prettyprinted on-screen representation of the polynomial. But, as long as you don't try to display it, Maple's kernel component (aka the "mathematical engine") performs operations on thousand-term polynomials nearly instantaneously and million-term polynomials in seconds. The following session transcript confirms this:

Download BigPoly.mw

So, what you should do is copy-and-paste your Mathematica polynomials to a plain-text file, for example a Notepad file if you're using Windows. I'm guessing that Mathematica's plain-text syntax is acceptable, but I'm not sure that it is.^[*1] Make the first characters of the file P:=  (or any other variable name). Make the last character of the file a colon (:). Then save that file, and use Maple's read command as shown above.

[*1] Primer on Maple's plain-text (generalized) polynomial syntax:

1. There is no implied multiplication via juxtaposition. All multiplication should use *. Do not use dot (.) for multiplication. Although it is acceptable to Maple syntax, it changes the meaning.
2. Constructions such as 2(x + y) are acceptable to Maple syntax, but they don't mean 2*(x+y)!
3. Exponentiation can be done with ^ or **. Although these mean exactly the same thing in Maple syntax, there can be no white space (white space includes line breaks) between the characters of a two-character operator such as **.
4. The only acceptable algebraic grouping symbol is round parentheses ( ). Mathematica's square brackets [ ] and curly braces { } are not correct. Although these will be accepted by Maple's syntax, they will change the meaning of the expression, and it won't be a polynomial.
5. Addition is +.
6. Integers are sequences of digits optionally preceded by + or -, but two operators cannot be juxtaposed (even with white space). So, for example x^(-1) is acceptable, but x^-1 and x^ -1 are not.
7. If you must include white space where it isn't otherwise allowed, precede each such character with backslash \. So, this can precede a line break.
8. Division is /.
9. Precedence of operators is exponentiation, followed by multiplication and division, followed by addition and subtraction.
10. *, /, +, - are associative left-to-right. ^ is not associative. So, for example, x/2/3 equals x/6, (x^3)^2 equals x^6, and x^3^2 is a syntax error.
11. (Ultimately, anything can be used as a variable name in Maple if it's properly quoted. I won't cover those rules here.) Any sequence beginning with a alphabetic character and followed by alphanumeric characters can be used as a variable. A few of these that have pre-defined meanings that you should be aware of are D, I, O, Pi, Beta, gamma, GAMMA, Beta, Zeta, and Psi.

1. Use Pi instead 3.14159 for the numeric solution. You should never use an explicit numeric approximation for Pi in Maple (unless you're writing an article about approximations of Pi).

2. Remove the inequalities from the analytic pdsolve command.

3. Define Ana via Ana:= unapply(rhs(sol), (x,t));

4. Define Sub via Sub:= (X,T)-> evalf(Ana(X,T) - eval(u(x,t), RT(X,T))); 

All of the techniques that you show (genfunc, rsolve) are for finding an exact recurrence pattern. The likelihood of being able to do that with data from nature is pretty slim. Instead, what you want to find is a statistical pattern. This a where the TimeSeries package might help. You should also change 1 (negative) to -1, 2 (positive) to 1, and 3 (neutral) to 0. "Neutral" not being between the poles is ridiculous.

Several points:

1. If you want standard matrix arithmetic, then you must declare them with Matrix (or its equivalent with rtable, etc.), not Array. The dot operator performs elementwise multiplication on Arrays.

2. Yes, for compiled code, you'd need to write explicitly the multiplications loops.

3. The arithmetic of dense matrices of hardware floats is already highly optimized in Maple. You'll not gain anything by compilation unless you can exploit some special structure particular to your matrices. For example, it's possible to speed up some operations on tri-diagonal matrices. (I wrote an Answer for this case (compiled code for tri-diagonal matrices). I put a link in a Reply.) 

4. Here are three ways to transpose a matrix or vector A:

1. LinearAlgebra:-Transpose(A)
2. A^+
3. A^%T 

There are many different ways that you could animate that. Here's one:

plots:-animate(
    plots:-shadebetween, 
    [sin(x+t), cos(x-t), x = 0 .. 2*Pi, positiveonly], 
    t= -Pi..Pi, frames= 200, gridlines, thickness= 3, labels= ["", ""]
);

Is that what you wanted?

In addition to what's already been mentioned by Tom Leslie and Thomas Richard, there are a few other potential problems with your ODE system. But if you can get the plot that you show, then you must've already overcome those problems. So, you can integrate p(x) as follows: You include one more simple ODE and one more initial/boundary condition directly in your system: 

solnum:= dsolve({ode, ics, diff(P(x),x) = p(x), P(0) = 0}, numeric);

Then, solnum(1) will show the values of the components of the system at x = 1. This incudes P(1), which is the integral that you want.

I recommend against including the options method and maxmesh in your dsolve command, unless you've already had an error message for the same system that suggests those options.

I see that after running dsolve, you've tried to use p as if it were a numeric function. This is possible, but the syntax required for it is slightly more complicated. If you still need to access that numeric function, let me know. But there's no need for this for plotting or integrating.

Virtually all computer languages, including Maple, allow for variable names that are longer than one letter. So xyz is simply a variable distinct from x, y, or z. In other words, Maple does not infer any mathematical relationship between xyz and x, y, or z.

It is sometimes possible in Maple to specify multiplication by putting a space between variable names, as in x y z. However, I don't recommend using this feature.

The Iterator package allows for truly sequential (as in one at a time) generation of combinatorial objects, including permutations. It is very efficient both in memory and time. For example, to generate all permutations of [10,J,Q,K,A] in a loop, you could do

for p in Iterator:-Permute(5) do 
    printf("%a\n", [10,J,Q,K,A][[seq(p)]]) 
od:

In this loop, the ps are generated as they are needed on each iteration rather than being generated as a whole by the Iterator:-Permute command itself. This is why the package is very memory efficient.

One aspect of the package that takes a little getting used to is that the objects are always returned as Vectors rather than lists or sets, which is why my index above is [seq(p)] rather than simply p. Even with the necessary conversions, the package is fast.

Do this:

`print/laplace`:= '`print/laplace`':
alias(X(s)= inttrans:-laplace(x(t), t, s)):

The first command gets rid of the fancy L (by dereferencing the procedure that changes the printing from plain laplace). The second command is a correction of your original alias command. You had the t and s reversed.

Okay, here is a procedure to do it (to solve your linear system). The results are returned as an n x s matrix each column of which is a vector akin to a solution vector of the ODE system. The vector b does not participate in this. The coefficient matrix from the ODE system can be constant (as in the example below from your PDF), or it can depend on the independent variable, t.

Num211:= proc(
   Mt::Matrix(square), gt::Vector[column], t::name, x0::Vector[column], 
   A::Matrix(square),
   c::Vector, h::numeric
)
uses LA= LinearAlgebra, AT= ArrayTools;
local 
    n:= upperbound(Mt,1), s:= upperbound(A,1),
    k, x, X:= Matrix((n,s), symbol= x),
    g:= unapply(gt,t), M:= unapply(Mt,t),
    f:= (t,x)-> M(t).x + g(t)
;
    ArrayTools:-Reshape(
        LA:-LinearSolve(
            LA:-GenerateMatrix(
                [seq(
                    AT:-Alias(X, 0, [n*s], 'column') 
                    - h*LA:-KroneckerProduct(A, LA:-IdentityMatrix(n))
                       . AT:-Concatenate(1, seq(f(h*c[k], X[..,k]), k= 1..s))
                    =~ AT:-Concatenate(1, x0$s)
                )],
                [seq(X)]
            )
        ),
        [n,s]
    )
end proc
:        
#Your example:
Num211(
    <1, 4; 3, 0>,  #M(t)
    exp(-2*t)*<4,3>,  #g(t)
    t,  #t
    <3,5>,  #x0
    <0, 0, 0; 1/4, 1/4, 0; 0, 1, 0>,  #A
    <0, 1/2, 1>,  #c
    0.00001   #h
);

            [3.  3.00013500083751  3.00027000335005]
            [                                      ]
            [5.  5.00006000093751  5.00012000375002]

This is a quick first draft. For large-scale usage, it could be made more efficient.

If the system of linear equations that you present has a unique solution (which is likely), then the final solution depends entirely on the s-vector b. You could choose b to give the exact solution for x(h), for example! In other words, figuring out a good b to use must be just as difficult as approximating the solution in the first place. Note that the system of equations does not involve b. Or, in more other words, the solution of the system of linear equations doesn't matter. So, either it's a crack-pot method, or your PDF is missing some crucial elements.

What you want cannot be easily done with implicitplot, but it can be easily done by parameterizing the curves and using animatecurve:

common:= color= red, thickness= 5, scaling= constrained, axes= none, paraminfo= false:
plots:-animatecurve([[-cos(t), sin(t), t= 0..Pi], [-cos(t), -sin(t), t= 0..Pi]], common);
plots:-animatecurve([[-3*abs(sin(y)), y, y= Pi..0], [-3*abs(sin(y)), y, y= Pi..2*Pi]], common);
 

The second solve command is missing a comma between its last two equations.

Note the error message:

Error, missing operator or `;`

This indicates a purely syntactic error rather than a mathematical error.

Your IC and BC are already sets, so {IC,BC}, which you pass to pdsolve, is a set of sets rather than a set of equations.

How about this?

restart:
interface(imaginaryunit= _i):
f:= (x,y)-> Re(sqrt(x+_i*y)):
plots:-display(
    plot3d(
        (common:= (
           [seq([x, y, k*f(x,y)], k= [-1,1])], x= -1..1, y= -1..1
        )), 
        grid= [200$2], style= surface, glossiness= 1, transparency= .05, 
        colorscheme= [
            "xyzcoloring", 
            [(x,y,z)-> y^2, (x,y,z)-> abs(x*y), (x,y,z)-> x^2]
        ]
    ),
    plot3d(
        common,
        grid= [25$2], style= wireframe, thickness= 3, 
        colorscheme= [
            "xyzcoloring",
            [(x,y,z)-> 1-y^2, (x,y,z)-> abs(x*y), (x,y,z)-> 1-x^2]
        ]
    ),
    lightmodel= light2, orientation= [55,75,0], projection= .85, axes= frame,
    labels= ['x', 'y', Re(w)]
);