Here is how:

See Maple's help on fieldplot3d.  This is an example from there:

with(plots):
fieldplot3d([2*x,2*y,1],x=-1..1,y=-1..1,z=-1..1,grid=[5,5,5]);

Your argument here is backwards: You should check the correctness of Matlab with the help of Maple, not the other way around.

Matlab does its calculations in the machine hardware which carries about 16 digits of precision.  Maple can do its floating point calculations in software with any degree of precision.  The default is 10 digits which is less than Matlab's, but you may request higher.  For instance:

K := 25888158.4218766;
K-sqrt(K^2-1);              # default 10 digits precision
    0.
evalf(K-sqrt(K^2-1), 30);   # 30 digits precision
    1.93138496702600*10^(-8)

Compare this correct answer to Matlab's not-so-good answer.

You may change Maple's default precision of 10 by setting the value of Digits , as in:
Digits := 30;
Then all floating point calculations will be done with 30 digists of accuracy in software.  This will slow things down, therefore I don't recommend doing that unless you really have a need for it.

As an alternative to what kitonum has suggested, you may try:

ss1:=-2:  ss2:=1/2:  ss3:=-5:
printf("The 3 equation in A, B, C for s=%a, s=%a, and s=%a:\n", ss1, ss2, ss3);
     The 3 equation in A, B, C for s=-2, s=1/2, and s=-5:

Here is an alternative to what rlopez has suggested.

restart;
with(Physics):
Typesetting:-Settings(typesetdot=true):

Here I will pick a Lagrangian function just for illustration.  Replace it with yours.

L := 1/2*a[1]*diff(q[1](t),t)^2 + 1/2*a[2]*diff(q[2](t),t)^2
   + b*diff(q[1](t),t)*diff(q[2](t),t)
   + c[1]*q[1](t) + c[2]*q[2](t);

The two Euler-Lagrange equations are:

de[1] := diff(diff(L, diff(q[1](t),t)), t) - diff(L, q[1](t)) = Q[1](t);
de[2] := diff(diff(L, diff(q[2](t),t)), t) - diff(L, q[2](t)) = Q[2](t);

Here is an alternative to the solution presented by vv and kitonum.  It involves a little bit of a "cheat" in embedding the drawings in 3D.  The method can handle non-convex domains with equal ease.  Here are a couple of samples extracted from the worksheet mw.mw.

You are attempting to solve your equations in two phases: first for p and q, and then for M.  It is easier to solve all of them together at once.

There are four unknowns, p1, p2, q1, q2, from the first set of equations.  There are 16 knowns, M[i,j], from the second set of equations.  Put them together and solve a system of 20 equations in 20 unknowns, and then you are done.

You didn't say what your differential equations are, therefore I will make up some arbitrary ones here.
de1 := diff(x(t),t,t) = -x(t) - diff(x(t),t) + y(t);
de2 := diff(y(t),t,t) = -y(t) - diff(y(t),t) + x(t);
ic := x(0)=0, D(x)(0)=1, y(0)=0, D(y)(0)=1;

As you see, I have isolated the second derivative terms on the left-hand sides of de1 and de2.

Now solve the system:
dsol := dsolve({de1, de2, ic}, numeric);

Let's say we want to examine the solution at time t=1.23.  We do
dsol(1.23);

The accelerations at time t=1.23 are obtained from
eval(rhs(de1), dsol(1.23));
eval(rhs(de2), dsol(1.23));

This is only a slight variation of a question that you asked a month ago.  Here is how to solve it:

Worksheet: mw.mw

restart;
P := <1,1,1> + s*<2,-1,-1> + t*<-1,3,2>;
Q := <1, -3, 5>/sqrt(35);
t_range := 0..1;
s_range := 0..1;
center := eval(P, {t=`+`(op(t_range))/2, s=`+`(op(s_range))/2});
plots:-display([
  plot3d(P, t=t_range, s=s_range, color=blue, style=surface),
  plots:-arrow(center, Q, color=red)
]);

_U1 is a "dummy" variable of integration, because the solution of your differential equation is expressed as a convolution.  You may replace _U1 by any other symbol you like, for example:
restart;
de := diff(x(t),t) = a*x(t) + b(t)*u(t);
dsolve({de, x(0)=alpha}, x(t), method=laplace);
subs(_U1=s, %);
                       

There may be shorter ways of doing this but this is what first came to my mind:

Download worksheet: mw.mw

In the attached worksheet I calculate the fastest path of a point mass sliding without friction on an arbitrary surface z=f(x,y) under the action of gravity which pulls in the negative z direction,. The path's endpoints are prescribed,

For simplicity I assume that the path can be expressed as a function of the x coordinate.  That's  the reason for the "partial" in the title.  The general case requires a fully parametric representation of the path.

Here is an illustration when the surface is z=f(x,y)=x²+y².

In the worksheet there is a formula for the overall travel time, however I don't know how to evaluate it because it is expressed in terms of an integral of an expression involving the solution of dsolve/numeric.  You guys may know how to do that.

brachistochrone-generalized.mw

Note added later

Following Carl Love's suggestion, I have added a section that calculates the travel time.  Here is the extended worksheet.

brachistochrone-generalized-ver2.mw.

The travel time for the illustrative example above is T=4.4019.