See ?printlevel . The default value of the environment variable printlevel is 1. Each level of loop nesting (or if statement nesting) adds one level. So, if printlevel is set to 2 or more, then you will see what's happening inside your double loop.

Local variables and parameters evaluate differently from global variables. Quoting from ?eval (seventh paragraph of Description):

 The default evaluation rules in Maple are full evaluation for global variables, and one-level evaluation for local variables and parameters.

When your eq2 (in procedure f) is evaluated to one level, that leaves its x unevaluated. So, the cure is to use eval(eq2), which forces a full recursive evaluation.

restart;
f:= proc()
local
     x,y,
     eq1:= 5+3*x=0,
     eq2:= 2+7*x-3*y-5*x*y=0
;
     x:= solve(eq1,x):
     y:= solve(eval(eq2),y):
lprint('x'=x,'y'=y);
end proc:

f();
x = -5/3, y = 29/16

To address your second question, about ||: The || operator simply was not intended to concatenate arbitrary expressions. It was intended as a convenient way to create indexed symbols. What you are trying to achieve can be easily done with the formatted print statement, printf:

printf("x is %a", -5/3);

The zip command is used to apply a two-parameter procedure to two matched arrays (or lists) of arguments.

restart:
B:= (x,y)-> fsolve(y = (1-exp(-x*b))/(1-exp(-50*b)), b):
X:= #Some array of x values
Y:= #Some array of corresponding y values
zip(B, X, Y);

Bisection is a method of last resort. Why do you want to use it? You can quickly get any number of roots of f(p) with RootFinding:-NextZero. After making Preben's unapply correction, do this:

N:= 100:  #Number of roots desired.
Root:= Array(1..N):
Root[1]:= -RootFinding:-NextZero(p-> f(-p), 0):
for n to N-1 do
     Root[n+1]:= -RootFinding:-NextZero(p-> f(-p), -Root[n])
end do:

The minus signs are because it seems that you want the roots in the negative direction.

All that you need is

series(y(x+h) - y(x-h), h=0, 3);

NumberToList:= proc(x::realcons)
uses S= StringTools;
     parse~(S:-Explode(S:-Select(S:-IsDigit, sprintf("%a", x))))
end proc:

Change the last line of your last for loop from

theta[ii](tau):=rhs(%);

to

theta[ii](tau):= value(rhs(%));

The value command makes Maple perform the integrations.

Since the desired simplification is not generally true for all bases, you have to use the symbolic option to simplify:

simplify(%, symbolic);

Using Maple 17, when I use solve followed by evalf, I get

Try doing a restart, then try the solve and evalf again.

A real answer can be obtained with fsolve or RootFinding:-NextZero in addition to the method that you mentioned.

KernelDensityPlot is a command in the Statistics package that will output the same plot in the normal inline manner:

Statistics:-KernelDensityPlot(Cdata);

AB,C:= LinearAlgebra:-GenerateMatrix(
     [F[k] $ k= 0..3],
     [u[k,n] $ k= 0..3, u[k,n-1] $ k= 0..3],
     augmented= false
);
A:= AB[.., 1..4];  B:= AB[.., 5..8];

1. In most cases, you reduce the error by increasing the value of Digits, which has the default value 10.

Digits:= 15:
fsolve(eqs2);

2. To get an exact solution to a decimal problem (if it's possible), use convert(..., rational):

solve(convert(eqs2, rational));

The midpoint method is not implemented in dsolve(..., numeric, method= classical[...]), so you'll have to write your own. Here's my ad hoc version---tailored to this specific problem. As you can see, the error plot is a straight (in log-log mode) line with a slope of nearly exactly 2.

Download Midpoint_error.mw

You need to change cos(kz - `ωt`) to cos(k*z - omega*t). I ran the code, and everything works if you correct this expression in the three places that it occurs.

I am not sure what you mean by "the scale of plot". Do you want a logarithmic scale? Do you want to extend the range on the y-axes so that both plots have the same y-axis?

Your "final result" has three terms with dimensions of force and one term with dimensions of force/temperature. Was there an integral in your code that you expected would cancel the temperature, which is your variables alfa1 and alfa2? If so, I couldn't find it.