The error message says it. The container objects in question are:
1. The matrix A having 4 elements.
2. The list LM having 2 elements.

However, this doesn't really explain the failure, because if LM contains 4 matrices you get the same error. Hopefully somewhere it is stated that the objects must be of the same type (in most cases?):
This doesn't work:
[1,2,3,4]*~A;
This does:
[1,2,3,4]*~[7,9,13,17]

@farazhedayati I believe I got it. See links above.

@farazhedayati Unfortunately, it appears that your two worksheets contain no code used for producing the graphs.
I suppose that you removed the code for presentation purposes.
I doubt that the graphs can be rescued from those worksheets. I hope I'm wrong.
##
I manage to get this:  contourplot.mw by some mouse work.
## One more: contourplot2.mw
## In the process I changed proportions, but you can adjust those by using the mouse.

Maybe someone can help you if you upload the worksheet.
The text looks like plot data. How it got there will (could) be revealed if we see your worksheet.

If your make a search in Maple for userinfo in the result from
showstat(limit);
you won't find userinfo.
Thus infolevel[limit]:=5;
wouldn't help anything.
You could try debug(limit); but that won't be very illuminating.
Or you could try stopat(limit); 
Basically, don't expect a computer program like Maple to work like a university teacher presenting things in a string of logically ordered steps as in a Calculus class.

Although the title could be seen as a question, it appears that it is not. The code seems to be running.
You ought to avoid implicit declarations of locals though.
Either declare those local or suppress the warning (which may not be that great an idea in general):
interface(warnlevel=1);

I don't know, but last year it was at the end of March.

To go from the inert Int to the active int, simply use value(res), where res is the result you got from dsolve.

The standard way to turn exact (symbolic) expressions like sqrt(2), Pi, sin(78), or what have you, into floating point numbers, is to use evalf:
evalf(Pi);
evalf(sin(78));
If you want to work with more digits than the default 10 (e.g. 15) you can do
evalf[15]( Pi );

@Markiyan Hirnyk An nth order ode needs n boundary (or initial) conditions.
If you try the following after your dsolve command you will see that ((D@@4)(w))(0) appears, i.e. an arbitrary constant:
indets(%,function);
### But using torabi's worksheet and adding method = laplace, as you did, works.

You must have introduced a floating point number (i.e. a decimal number) at some stage. After that things like you describe can happen.
To point out where it happened we need to see the code. You could upload a worksheet.

The site was down for maintenance a day or so ago, but it surely didn't last for days.

Do you need to refresh the page or clear the cache in your browser?

I have been able to access the page each day for months with the exception mentioned. Needless to say, I'm not looking at it at all times though.

@kuker2008 You could use infnorm from the numapprox package. It finds the infinity norm of the expression, i.e. the supremum (or maximum) of the absolute value of the input on the given interval.
Thus there is no difference between these two:
 

numapprox:-infnorm(abs(b-a),x=-1..1);
numapprox:-infnorm(b-a,x=-1..1);

@Markiyan Hirnyk So Markiyan, if that is right, then how do you explain the difference between the very different results of

int(int(F,x=-infinity..infinity),y=-infinity..infinity);

and

int(F,[x=-infinity..infinity,y=-infinity..infinity]); 

where F is given by

f:=x->1/sqrt(2*Pi)*exp(-x^2/2);
F:=convert(f(x)*f(y)*x*x*abs(x+y),piecewise,x);

And is the note I quote from the help page for int in my answer false in its statement about a difference?

##### But I notice that my argument for VectorCalculus:-int using the iterated form was faulty.
Consider this:
 

restart;
int(phi(x,y),[x=-infinity..infinity,y=-infinity..infinity]);
lprint(%);

The result of lprint is the iterated form. I take that as an indication that after the "assumptions" job is done the integral is computed via an iterated integral.
That is in full agreement with the help page note.

First of all you need to use a midpoint method as the error message says, when you execute the dsolve command.
If you solve for the highest derivatives there is a denominator in the ode for H of the form eta*(....). So that will cause problems at eta=0.
Thus change that dsolve command to
dsys1 := {Eq1, Eq2, Eq3, Eq4, IC1, IC2}; dsol1 := dsolve(dsys1, numeric, output = listprocedure, method=bvp[midrich]);
or use the other midpoint method bvp[middefer].
I tried with no luck whatsoever.
Do you have any reason to believe that your system has a solution at all?