Let's note that 

u := (1 + 20230321)*x*y - (x^2 + y^2)/2

is an homogenous polynomial wrt x and y.
One can rewrite u by setting [x=lambda*cos(t), y=lambda*sin(t)] which gives 

v := collect(eval(u, [x=lambda*cos(t), y=lambda*sin(t)]), lambda);
    /                         1       2   1       2\       2
    |20230322 cos(t) sin(t) - - cos(t)  - - sin(t) | lambda 
    \                         2           2        /        

The maximum (maxima in reality) of v are obtained when lambda is equal to +/- infinity and when the t unction (op(1, v)) reaches its maxima:

v       := collect(eval(u, [x=lambda*cos(t), y=lambda*sin(t)]), lambda):
opt_t_1 := Optimization:-NLPSolve(op(1, v), t=0..Pi, maximize):
opt_t_2 := Optimization:-NLPSolve(op(1, v), t=Pi..2*Pi, maximize):

opt_x_1 := op(2, v) *~ eval([cos(t), sin(t)], opt_t_1[2]):
opt_x_2 := op(2, v) *~ eval([cos(t), sin(t)], opt_t_2[2]):

identify(evalf[8](opt_x_1));
identify(evalf[8](opt_x_2));

              [1  (1/2)       2  1  (1/2)       2]
              [- 2      lambda , - 2      lambda ]
              [2                 2               ]

            [  1  (1/2)       2    1  (1/2)       2]
            [- - 2      lambda , - - 2      lambda ]
            [  2                   2               ]

Then the maxima are rejected to x = +/- infinity and y = x.

Which s a far better result than Mathematica provides.

So, let me ask you this question (and @QM too): suppose your are part of a MMA users blog and MMA tells you that 

 {\[Infinity], {x -> Indeterminate, y -> Indeterminate}}

wouldn't you say that MMA doesn't give a completely satisfactory answer ?

I feel a sense of unfairness in the question and this two-man bashing of Maple, am I wrong? 

d1 is a centered divided difference approximation of diff(u, x$3) at "point" [i+1/2, j];not an approximation of the laplacian of u at [i, j]
 

Download DividedDifferences.mw

This code allows you to print your list in a way that is close to what you want.
Printed lines have arround n=10 characters.

Of course arrangements have to be done if some elements of the list are numbers with more than n digits.
More of that, the list is converted into a string before being parsed, which means that spaces and commas are included in the calculation of the string to plot.
But I think this can give you some hints to go further (for instance if you want to plot numbers, you must augment the the range to plot 1..pc[cut] by 2*p, where p is number of occurrences of a comma in newLs[1..pc[cut]].

 

Download SmartPrint.mw

Unless I'm mistaken you miss one equation:

• eq1 is of the form

   A*diff(h(r), r$4) + B*diff(h(r), r$4) + LowerOderTerms = 0

  A and B do not contain derivatives of h anq q of order 4

• eq2 is of the form

   C*diff(h(r), r$3) + E*diff(h(r), r$3) + LowerOderTerms = 0

  C and E do not contain derivatives of h anq q of order larger than 2

• Then plugging eq2 into eq1 produces a single ODE of the formeq1 is of the form

   P*diff(h(r), r$4) + Q*diff(h(r), r$4) + LowerOderTerms = 0

  where P and Q do not contain derivatives of h anq q of order larger than 3

So you miss one ODE.
The detailed explanation is given in the attached file.
A toy problem is also provided to help you understand why dsolve doesn't work asyou expected.
 

Download dsolve_mmcdara.mw

A toy_problem dsolve can't find the solution of, but whose solution can be obtained after reparameterisation:

 

Download Toy_Problem.mw

restart:
sys := {x^2*y+x=1, x*y-y^2=0}:
solve(sys);
    {x = 1, y = 0}, 

       /          /  3         \            /  3         \\ 
      { x = RootOf\_Z  + _Z - 1/, y = RootOf\_Z  + _Z - 1/ }
       \                                                  / 
infolevel[solve]:=100:
solve(sys)
Main: Entering solver with 2 equations in 2 variables
Main: attempting to solve as a linear system
Main: attempting to solve as a polynomial system
Main: Polynomial solver successful. Exiting solver returning 1 solution
    {x = 1, y = 0}, 

       /          /  3         \            /  3         \\ 
      { x = RootOf\_Z  + _Z - 1/, y = RootOf\_Z  + _Z - 1/ }
       \                                                  / 

 

This suggests solve calls some polynomial solver.
Searching for it in the help pages returns SolveTools/PolynomialSystem. I don't know if solve calls this procedure (showstat(solve) is hard to browse). Assuming it does:
 

SolveTools:-PolynomialSystem(sys, {x,y});

Main: polynomial system split into 1 parts under preprocessing
Main: trying resultant methods
PseudoResultant: normalize equations, length= 43
PseudoResultant: 353502 [1 3999999986 y 1] 2 2 43 0 3 0
PseudoResultant: normalize equations, length= 13
PseudoResultant: normalize equations, length= 35
PseudoResultant: 30000 [1 100000559 x] 1 1 13 0 3 1
PseudoResultant: normalize equations, length= 3
PseudoResultant: -10 [] 0 0 3 0 3 1
PseudoResultant: 201729 [1 1000012108 y] 2 2 35 1 8 0
PseudoResultant: normalize equations, length= 21
PseudoResultant: 117197 [1 700003740 x] 1 1 21 1 8 0
PseudoResultant: normalize equations, length= 3
PseudoResultant: -10 [] 0 0 3 0 3 0
PseudoResultant: 2 solutions found, now doing backsubstitution
PseudoResultant: backsubstitution of x
PseudoResultant: backsubstitution of y
PseudoResultant: backsubstitution of x
PseudoResultant: backsubstitution of y
    {x = 1, y = 0}, 

       /          /  3         \            /  3         \\ 
      { x = RootOf\_Z  + _Z - 1/, y = RootOf\_Z  + _Z - 1/ }
       \                                                  / 

restart
p := 35*x^8 + 80*x^7 + 140*x^6 + 504*x^5 + 1260*x^4 + 1960*x^3 + 1960*x^2 + 1120*x + 10:
dp := diff(p, x):
fsolve(dp);
                         -1., -1., -1.

But I suspect that you are looking for an answer you can obtain "by hand" aren't you?
Here is, written in Maple, a step-by-step procedure which only requires doing polynomial division (which can be done by hand)
Download ByHand.mw

Details:

• Observe -1 is a root of dp = diff(p, x)
• Thus dp = Pol(x) * (x+1)^r with r >=1 and degree(Pol) <= 6
  Find the value of r this way:
  1. Set Q =Pol(x)
  2.  Divide (can be done by hand) Q by (x+1) and note R the remainder
  3. Go to point 1 until R = 0
• The value of r is the number of times you executed the steps 1 to 3
  Here r=3 ==> dp = Pol(x) * (x+1)^3 with degree(Pol) = 4
• Search if Pol(x) has a real zero Z by applying the procedure given at points 1 to 3.
  Here you recursively divide 4 times Pol(x) by (x-Z) : if the last remainder is 0, then is a root of dp, if not (which is the case) dp has only 1 real root (-1) of multipliciy 3.

Thus p has a only one stationnary point at x=-1 (I don't know if you must answer 1 or 3 ???)

The goal is to produce almost exactly the expressio of the solution you give, thus we have to go beyond the dsolve to reshape the solution Maple provides.
The procedure is described step-by-step

Download OneWay.mw

(1/2, 1/2) is not always a solution, look here

Download solution.mw

Maybe you missed something?

Probably a hideous way to proceed

Using  result := ssystem("tasklist <opts>"), where < opts >depends on what OS you use (Windows only) you can get a table (iin Windows sense) which contains all the PIDs.
This require to "postprocess" resultto display only the PIDs associated to a mserver.exeprocess

When you run xmaple, open a worksheet run the previous command to get the mserver's PIDs.
Proceed this way ach time you open a new worksheet: by comparing the new PID's list to the one of the previously opened worksheet you will get the PID associated to your current worksheet.
This can be useful is the stopbutton doesn't cancel the computations and if you can't close the current document.

Unfortunately I don't know how to identify to what node a mserverprocess belongs to when I distribute the computations over several nodes while using the GridPackage (being able to kill the computations on a specific node is something I would really like to have, all the more so that even after leaving Maple the process still exists and keeps using memory).

The command which generates output (4) is incorrect:

• You use = instead of := 
• You use exp^(...) instead of exp(...) 

Here is the correction (done with Maple 2015 whichh doesn't accept this awfull notation C(x, t) := ...)

 

Download easy_way_mmcdara.mw

L:=[[3, 2], [2, 1], [1, 2], [1, 2], [2, 3], [2, 1], [1, 2], [1, 1], [2, 1], [1, 2], [1, 1], [2, 1], [1, 1], [1, 3], [1, 2], [2, 1], [1, 3], [1, 3], [1, 3], [1, 2], [2, 2], [2, 3]]: 

Ls := convert~(L0, set): 
L1 := map(u -> [[u[]], ListTools:-Occurrences(u, Ls)], {Ls[]}): 
Lk := map(u -> if numelems(u[1])=1 then [[u[1][1]$2], u[2]] else u end if, L1)

{[[1, 1], 3], [[1, 2], 11], [[1, 3], 4], [[2, 2], 1], [[2, 3], 3]}

L2 := [map(u -> [add(op(u[1])), u[2]], Lk)[]]

           [[2, 3], [3, 11], [4, 1], [4, 4], [5, 3]]

answer :=(`*`@op)( (`*`@op)~(L2))

                             190080

If you want to find the argument which minimizes the sum of the norms for a given value of lambda, this is quite simple:

w[opt] := (X^+ . X + lambda*~J)^(-1) . X^+ . Y

where X is a N by P matrix, Y a column vector of length N and J the identity matrix of dimension N.

But maybe you are talking about Ridge regression where you seek fot the couple (w[opt], lambda[opt]) which minimizes thz sum of the norms???

Several things are weird or incorrect:

• Do you really need this weird Digits:=100???
  It works quite the same with Digits:=10.
   
• p(t) is a nonsense, for the integration is done over t and thus its result cannot depend on t:

  p := int(f(t)*exp(-lambda^2*t), t=0..L)
                                       -134
                        -4.804438634 10    
  

• As p is a constant, C(x, t) is a function of x alone.
  The message 
  Error, (in sum) summation variable previously assigned, second argument evaluates to 200. = 1 .. 500
  is clear: you have previously assigned n to 200.
  The cure: protect n by writing 'n' =1..500

  C := x -> sum((2/L*sin(n*Pi*x/L)*exp(1)*sin(n*Pi*x/L))*p, 'n' = 1 .. 500):
  C(1);
                                       -136
                        -1.871559740 10    
  

With these correction the worksheet doesn't contain any error...
But are you sure that the integral

int(f(t)*exp(-lambda^2*t), t=0..L)

is the one you wanted to write???

In Maple 2015.2 SumTools doesn't offer a procedure Change as IntegrationTools does.
So the idea to transform the sum/Sum into Int, to apply Change, and to go back to the sum/Sum.

S := sum(a[k]*(k+r)*(k+r-1)*x^(k+r-1), k = 0 .. infinity): 
J := Int(op(S)): 
IntegrationTools:-Change(J, k=j+1, j):
IntegrationTools:-Change(%, j=k, k): 
Sum(op(%))
         infinity                                     
          -----                                       
           \                                          
            )                                  (k + r)
           /     a[k + 1] (k + 1 + r) (k + r) x       
          -----                                       
          k = -1                                      

@sursumCorda 

You're right.
Here is another way based on the following principle:

• Let A some term in rules.
• Set B = eval(A, R) where R denotes the rules with element removed
• Keep proceeding this way while B <> A after setting R=R \ B and A=B.

I initially removed elements of the form b=b (trivial 1-cycles).
 

Download cycles_1.mw


EDITED
to get all the cycles without permutations removed, replace the last for..do..end do structure by
 

# AEC stands for All Elementary Cycles (rotations removed) 

AEC := NULL:
for i from 1 to numelems(rules) do
  LocalCycle := cycle(i):
  if `not`(has(convert~([AEC], set), {convert(LocalCycle, set)})) then
    AEC := AEC, LocalCycle;   end if:   printf("from = %a : cycle = %a\n", rules[i], LocalCycle);
end do:
from = g = d : cycle = [g = d, d = c, c = g]
from = e = a : cycle = [e = a, a = e]
from = a = b : cycle = []
from = f = d : cycle = [f = d, d = c, c = f]
from = c = g : cycle = [c = g, g = d, d = c]
from = d = c : cycle = [d = c, c = g, g = d]
from = c = f : cycle = [c = f, f = d, d = c]
from = a = e : cycle = [a = e, e = a]
from = f = g : cycle = [f = g, g = d, d = c, c = f]

print~({AEC} minus {[]}):
                         [e = a, a = e]
                     [f = d, d = c, c = f]
                     [g = d, d = c, c = g]
                  [f = g, g = d, d = c, c = f]