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)]
);

Here's a small wrapper procedure for DEtools:-matrixDE that does everything that that does, then solves for the integration constants given the initial conditions, and then returns the solution in Matrix and Vector form. The wrapper accepts the options used by maxtrixDE---method and solution---if they are placed at the end of the command.

Download MatrixDEwithIC.mw

My procedure also converts the output from matrixDE from the cruddy old-style array format to the modern Matrix and Vector format.

Your expected answers are transposes of what they should be, right?

And you seem too have ignored the advice from the Answers to your previous Question.

Anyway, the transposes of your expected answers can be obtained like this:

restart:
n:= 2:
X:= Vector(n, symbol= x);
W:= Matrix(n$2, symbol= w);
V:= W.X;
seq(diff~(v,W), v= V);

Of course, these can be transposed if you want.

Continuing with Preben's fine advice, you should investigate initial conditions x(0) = c, where c is between a and b/2. This is the middle band, where the more-interesting (and usually more-practical) solutions lie. 

The issue is not that you're trying to differentiate a matrix. Rather, the issue is that you're differentiating with respect to x, a vector, not a single variable. Use Gradient instead of diff:

restart:
x:=<x1,x2>;
W:=<<w11, w12>|<w21, w22>>;
Cal:= (W.x)^%T.(W.x);
VectorCalculus:-Gradient(Cal, [seq(x)]);

Note that I've made x a vector rather than a matrix, and this makes Cal a scalar rather than a one-element matrix.

The issue is simply a misplaced parenthesis, as pointed out by Tom. The purpose of this Answer is to show you two ways to construct L more robustly so that perhaps you would've been able to figure that out just from the error message. Your L:= proc(N) N end is not really an identity function; rather, it's a projection function of its first argument. Note that L(1,2) returns 1. Here are two ways to make L better:

1. L:= proc(N, $) N end proc: This makes L a true one-argument identity function: Any attempt to pass L more (or less) than one argument will give an error message. So, the error message from your plot command would've been:

Error, invalid input: too many and/or wrong type of arguments passed to L; first unused argument is k = 0 .. 10

I think that it's obvious from that message that misplaced parentheses are the issue.

2. L:= ()-> args: This makes L a true identity function of an arbitrary number of arguments: All arguments are returned. The error message would've been:

Error, (in plot) incorrect first argument [min(5, max(-5, 5^(10*x)+5^(x*k))), k = 0 .. 10]

In this case, the special evaluation rules of seq prevent it from seeing the k= 0..10 as its optional second argument because that expression is part of the result from L. L hasn't yet been evaluated; that's the "special" part. So, seq sees the call to L as its sole argument, and it does what it usually does when it only has one argument, which is kind of trivial in this case: It returns the operands of L, which are the arguments. It's just a coincidence that this is what L would've done anyway. Those special evaluation rules could be made normal by replacing seq with (value@%seq). I'm not suggesting that you do this, because it's obvious that you didn't intend for k= 0..10 to be passed to L anyway. Nonetheless, this modification does give the correct plot.