Problems of this type can be done by finding the shortest path in a graph, although I don't know if this method has the best asymptotic efficiency. For this particular problem, the answer 11 (11 transitions) is returned instantaneously, alongs with the transitions themselves.

The vampires, maidens, and the lift are a system. We model the states of the system by who and what is in the lobby after the completion of each movement of the lift. (What's in the lobby necessarily determines what's in the bar.) Each state is represented by an ordered triple [V,M,L], where V is the number of vampires in the lobby, M is the number of maidens in the lobby, and L is the number of lifts in the lobby. We are going to model the system as a directed graph. Its vertices will be the valid states. Hence the name StatesVx. First, determine the valid states.

StatesVx:= table():
for V from 0 to 3 do
     for M from 0 to 3 do
          for L from 0 to 1 do
               if
                   #Can't have more V than M in lobby unless M=0:
                   (V > M implies M=0) and
                   #Can't have more V than M in bar unless all M are in lobby:
                   (M > V implies 3=M) and
                   #L can never be left where there are no people:
                   (V+M=0 implies L=0) and (V+M=6 implies L=1)
               then
                    #Maple won't allow lists as vertex labels, so I convert them to names:
                    StatesVx[V,M,L]:= nprintf("%a", [V,M,L])
               end if
          end do
     end do
end do;

States:= [indices(StatesVx)]:

The arcs of the graph will be the transitions of the states that can be achieved with a single movement of the lift. We examine each ordered pair of states and determine whether it's a possible transition.

Transitions:= table():
for pair in combinat:-permute(States,2) do
     Vx||(1..2):= pair[]; (V1,M1,L1):= Vx1[]; (V2,M2,L2):= Vx2[];
     if
          #The lift must move:
          L1 = 1-L2 and
          #If the lift comes down, the number of people in the lobby increases:
          (L1=0 implies V1+M1 < V2+M2) and
          #If the lift goes up, the number of people in the lobby decreases:
          (L1=1 implies V2+M2 < V1+M1) and
          #The next 4 conditions say that if V in the lobby increase,
          #then M can't decrease, and vice versa:
          (V1 > V2 implies M1 >= M2) and
          (V2 > V1 implies M2 >= M1) and
          (M1 > M2 implies V1 >= V2) and
          (M2 > M1 implies V2 >= V1) and
          #The total change of people in the lobby is at most 2:
          abs((V1+M1) - (V2+M2)) < 3
     then
          Transitions[StatesVx[Vx1[]],StatesVx[Vx2[]]]:= 1
     end if
end do:

Now compute a shortest path through the graph from the initial state (everything in the lobby) to the final state (nothing in the lobby).

GraphTheory:-ShortestPath(
     GraphTheory:-Graph({indices(Transitions)}),
     StatesVx[3,3,1], StatesVx[0,0,0]
);

     [`[3, 3, 1]`, `[1, 3, 0]`, `[2, 3, 1]`, `[0, 3, 0]`, `[1, 3, 1]`, `[1, 1, 0]`, `[2, 2, 1]`, `[2, 0, 0]`, `[3, 0, 1]`, `[1,        0, 0]`, `[1, 1, 1]`, `[0, 0, 0]`]

So the required number of transitions is one less than the number of vertices in that shortest path.

nops(%)-1;

     11

Download Vampires_&_Maidens.mw

Suggestions for improvement of the above algorithm are most welcome, especially improvement of the asymptotic efficiency as the number of vampires and maidens increase. Solutions that merely count the number of transitions without listing them are not welcome (as comments to this Answer).

You need to use a units environment. There are two pre-defined environments, Standard and Natural.

restart:
with(Units):
with(Natural):
UseSystem(SI);
phi:= 30*arcdeg;

sin(phi);

This is suprisingly simple in Maple:

diff(ln(x)/x, x$n);

diff(exp(x)*ln(x), x$n);

Yes it is easily done via

plot([<X|Y1>, <X|Y2>])

or, if you have a large numbers of Y's, you may want to map the adjoining of X, like this:

plot(map2(`<|>`, X, [Y||(1..2)]))

The following is a close match to the plot in your Question:

X:= <0,1,2,3>:
Y1:= <0.5, 2.5, 1.75, 2>:
Y2:= <0.3, 1.4, 1.6, 3.2>:

plot(
     map2(`<|>`, X, [Y||(1..2)]),
     color= ["Orange", "LightBlue"],
     thickness= 4, view= [-0.5..3.5, 0..3.5],
     axes= frame,
     axis[2]= [gridlines, tickmarks= [seq(.5*k, k= 0..7)]]
);

Maple's command for solving differential equations is dsolve, not solve. If Maple could solve your equation, the command to do so would be

dsolve(diff(y(x),x,x) = lambda*x*y(x)/sqrt(-1+x));

(Note that there're no square brackets and no capital letters in that. Note also how the second derivative is specified.) But Maple (Maple 16 at least) cannot solve this equation. The above command returns a DESol structure, which is just an abstract representation of an unknown function.

Use the command odetest.

First note that in differential equations, derivatives must be entered as diff(E(r),r), not diff(E,r). Here's your system of ODEs, corrected:

odesys:= {
   diff(E(r),r)^2+diff(B(r),r)/E(r)+E(r)*B(r)*r=0,
   diff(B(r),r)^2/E(r)+5*r = 0
}:

And here's your putative solution:

puta:= {E(r) = r^a, B(r) = exp(a*r)}:

odetest(puta, odesys);

The fact that the above is not identically (with respect to r) zero for any a shows that puta is not a true solution.

I only have time right now to give a brief example of how this might be accomplished.

The interpolation will be handled by command CurveFitting:-ArrayInterpolation (see ?CurveFitting,ArrayInterpolation). In the example, I construct arrays that will interpolate an approximation to cosine. Then I write a procedure that does the call to ArrayInterpolation. Then I use this procedure in a differential equation whose solution should be an approximation of sine. The interpolating procedure must be declared to dsolve as known.

Download known.mw

@Markiyan Hirnyk 

To avoid catching expressions with integer exponents, you need a more-specific type specification. This works for any expression ex:

indets(ex, algebraic^fraction)

The sqrt(-1), aka I, needs to be treated separately. If you want to catch this also, use

indets(ex, {identical(sqrt(-1)), algebraic^fraction}),

which will work even if the imaginary unit has been redefined.

Complicated roots are sometimes expressed in RootOf form. To catch these also, use

indets(ex, {specfunc(polynom, RootOf), identical(sqrt(-1)), algebraic^fraction}).

Here's something to try; I don't know if it'll work. Open the file in Maple 16, and then save the expression using the command

save expr "s.maple.m";

The .m causes the file to be saved in Maple internal format. This may be easier for Maple 18 to read.

What I see from a quick read through is that you never define r1 or r2. You can't plot them if you don't define them.

Here's a basic module factory to accumulate arithmetic means (not exponential means). I just wanted to show the basics of a module factory, and how simple they can be.

Download module_factory.mw

Answering your first question:

local Pi;
cos(x+Pi);

Regarding your second question: I don't know. The answer is not well defined. I mean, why should it return cos(x+Pi) rather then cos(x-Pi)? They're both equal to -cos(x).

There is no need for Arrays as long as you never assign to indexed list entries, which is not allowed for lists longer than 100 elements. By using seq to create lists and map to operate on them, you never need to assign to list entries. Here is a vastly simplified version of your code that produces exactly the same output and works for long lists.

Download plotScript.mw

Your integral equation can be converted to an ODE IVP. Then Maple's numerical ODE solver can be used on it.

Download Int_eqn_to_ode.mw

I am only debugging here the expression which updates z2 because the others were commented out. The errors in the commented-out expressions are very similar.

There are two bugs, both of them in the last factor in your expression for updating z2. First, you have arro[k] immediately in front of the parentheses that surround the (cosine + I*sine) expression. When a name is immediately followed by a left parenthesis, it is interpretted as function application rather than multiplication. You should put an explicit multiplication sign between the arro[k] and the left parenthesis. 

The second bug is that all function applications do require parentheses. That includes sin and cos, even though the parentheses are often omitted in standard mathematical notation and typesetting. So correcting both errors, the final expression may be entered as

arro[k]*(cos(2*Pi*j/21)-I*sin(2*Pi*j/21))