For "generic" values of beta, the only solution is the trivial one.  The eigenvalues are the values of beta^4 for which there is a nontrivial solution.  By scaling, the eigenfunction could be taken as one with D(y)(0) = 1.  Thus I would do it this way:

> myBCs:= y(0) = 0, D(y)(0) = 1, (D@@2)(y)(0)=0,y(L)=0;
  S:= dsolve({Y_diffEq, myBCs});
  condition := eval(diff(rhs(S),x,x),x=L) = 0;
  

Pi*_Z2/L*I, Pi*(2*_Z3+_B2)/L

i.e. the solutions are beta = n*Pi/L and beta = n*Pi*I/L for integers n (all of which have beta^4 = n^4*Pi^4/L^4).

In addition to Joe's excellent points: your Step3 ends in Step2, but according to your flow chart it should be followed by  UnitCk.

There are many mistakes here.

1) Your if and z -> are mixed up.  The first argument to selectremove is supposed to be a procedure.  That procedure can be of the form z -> ... .   The body of the procedure, the part after the ->, could be an if ... end if statement (except that, due to a bug, 2-D Math input has trouble with this).  But actually you don't need an if statement at all here.

2) The then clause in an if must come before the else clause.

3) The result of the procedure must be true,  false or fail.  Equations or inequalities can be evaluated to boolean with evalb or is (if non-numeric constants such as Pi or sqrt(2) could be involved, you want is).
 

4) Real part is Re, not R.

4) Your logic is wrong.  You said "R(z)>=0 or if R(z)=0 then Im(z)>0".  But if R(z) >= 0, the second part is not needed, and if R(z) < 0 the second part is true (because "A implies B" is true if A is false).  I suspect you mean "R(z) > 0 or (R(z) = 0 and Im(z) => 0).  Thus:

> selectremove(z -> is(Re(z) > 0 or (Re(z) = 0 and Im(z) >= 0)), %);

Note that the equations are symmetric in x1 and x2.

> Q:=eval(-3*x1 + 4*x1^3 - 3*x2 + 4*x2^3, x2=m-x1);
  plots[implicitplot](Q, x1 = 0..1, m = 0.4 .. 1);

This is actually a part of the ellipse 4 m^2 - 12 m x1 + 12 x1^2 = 3.

Note that it is only for sqrt(3)/2 <= m <= 1 that you get both x1 and x2 in the interval [0,1].

Getting the equations using "Image Properties":

>  e1:= 55200.0*(1+(3.20000000*10^7*a+1.200000*10^5*b+400*c+d)^2)^.5*
(2.048000000*10^15*a+7.680000000*10^12*b+2.560000000*10^10*c+6.40000000*10^7*d)
/((3.20000000*10^7*a+1.200000*10^5*b+400*c+d)*(4.800000*10^5*a+1200*b+2*c))
-1.177600000*10^12*(1+(3.20000000*10^7*a+1.200000*10^5*b+400*c+d)^2)^1.5
/((3.20000000*10^7*a+1.200000*10^5*b+400*c+d)^2*(4.800000*10^5*a+1200*b+2*c))
-1.766400000*10^10*(1+(3.20000000*10^7*a+1.200000*10^5*b+400*c+d)^2)^1.5
/((3.20000000*10^7*a+1.200000*10^5*b+400*c+d)*(4.800000*10^5*a+1200*b+2*c)^2) = 0; 
> e2:= 55200.0*(1+(3.20000000*10^7*a+1.200000*10^5*b+400*c+d)^2)^.5*
(7.680000000*10^12*a+2.880000000*10^10*b+9.60000000*10^7*c+2.400000*10^5*d)
/((3.20000000*10^7*a+1.200000*10^5*b+400*c+d)*(4.800000*10^5*a+1200*b+2*c))
-4.416000000*10^9*(1+(3.20000000*10^7*a+1.200000*10^5*b+400*c+d)^2)^1.5
/((3.20000000*10^7*a+1.200000*10^5*b+400*c+d)^2*(4.800000*10^5*a+1200*b+2*c))
-44160000*(1+(3.20000000*10^7*a+1.200000*10^5*b+400*c+d)^2)^1.5
/((3.20000000*10^7*a+1.200000*10^5*b+400*c+d)*(4.800000*10^5*a+1200*b+2*c)^2) = 0; 
> e3:= 55200.0*(1+(3.20000000*10^7*a+1.200000*10^5*b+400*c+d)^2)^.5
*(2.560000000*10^10*a+9.60000000*10^7*b+320000*c+800*d)/((3.20000000*10^7*a
+1.200000*10^5*b+400*c+d)*(4.800000*10^5*a+1200*b+2*c))-14720000*(1+
(3.20000000*10^7*a+1.200000*10^5*b+400*c+d)^2)^1.5/((3.20000000*10^7*a
+1.200000*10^5*b+400*c+d)^2*(4.800000*10^5*a+1200*b+2*c))-73600*(1+
(3.20000000*10^7*a+1.200000*10^5*b+400*c+d)^2)^1.5/((3.20000000*10^7*a
+1.200000*10^5*b+400*c+d)*(4.800000*10^5*a+1200*b+2*c)^2)+18400*(1+d^2)^1.5
/(c^2*d) = 0;
> e4:=  55200.0*(1+(3.20000000*10^7*a+1.200000*10^5*b+400*c+d)^2)^.5*
(6.40000000*10^7*a+2.400000*10^5*b+800*c+2*d)/((3.20000000*10^7*a
+1.200000*10^5*b+400*c+d)*(4.800000*10^5*a+1200*b+2*c))-36800*(1+
(3.20000000*10^7*a+1.200000*10^5*b+400*c+d)^2)^1.5/((3.20000000*10^7*a
+1.200000*10^5*b+400*c+d)^2*(4.800000*10^5*a+1200*b+2*c))-55200.0*(1+d^2)^.5
/c+18400*(1+d^2)^1.5/(c*d^2) = 0;
> eqs:= convert({e1,e2,e3,e4},rational);
> S:= solve(eqs);
  map(allvalues, [S]);

[{a = a, b = b, c = -80000*a-300*b, d = I}, {a = a, b = b, c = -80000*a-300*b, d = -I}, {a = a, b = b, c = -1/200*I-80000*a-300*b, d = I}, {a = a, b = b, c = 1/200*I-80000*a-300*b, d = -I}, {a = 1/32000000*2^(1/2)+1/80000*c+1/16000000*I, b = -1/150*c-1/80000*2^(1/2)-1/40000*I, c = c, d = 1/2*2^(1/2)}, {a = 1/32000000*2^(1/2)+1/80000*c-1/16000000*I, b = -1/150*c-1/80000*2^(1/2)+1/40000*I, c = c, d = 1/2*2^(1/2)}, {a = -1/32000000*2^(1/2)+1/80000*c+1/16000000*I, b = -1/150*c+1/80000*2^(1/2)-1/40000*I, c = c, d = -1/2*2^(1/2)}, {a = -1/32000000*2^(1/2)+1/80000*c-1/16000000*I, b = -1/150*c+1/80000*2^(1/2)+1/40000*I, c = c, d = -1/2*2^(1/2)}]

Actually the last four are not really solutions, because they make 480000.0000*a+1200*b+2*c = 0 and thus lead to a division by 0 in the equations.

Note that using EDP[2] your IBC[2] reduces to i(100)*(9*s(100)-4) = 0.  There are really two possibilities: i(100) = 0 or s(100) = 4/9.  Since your IBC[1] says s(100) = 1, you can't have s(100) = 4/9, so the second boundary  condition is really i(100) = 0.
And then dsolve(..., numeric) will happily solve your system, but of course the solution is s(z) = 1, i(z) = 0. 

No trick, just a different command: listdensityplot.

You might use the axis option.  For example, axis=[location=low] will make each axis's tickmarks be on the edge corresponding to the lowest values of the other two variables.

In general there's no way to do this (i.e. to retrieve the command that produced a certain result).  There is an Undo item in the Edit menu (with keyboard shortcut Ctrl-Z), but that is grayed out after you enter a command. You might also be able to use Restore Backup (File, Recent Documents, Restore Backup)  if Maple happened to auto-save the worksheet while the command was still there.

Enter (in a Standard worksheet): 

> 2^40000;

15842603725730786800597361511643477938585056174993108430408926417163696577122\
  88874451832215496349004[...11842 digits...]64589958897311815153669921216467\
  99484887367079871289697795984721749551119434590853334711885025509376

The result is 12042 digits long, but Maple shows only the first  100 and last 100, with [...11842 digits...] in the middle.  That's digit elision.  The settings determine the smallest number of digits that will trigger digit elision, and how many initial and final digits will be shown. Thus with the default settings 2^30000 (which has only 9031 digits) would be shown in its entirety, but if you set Threshold to 1000, Leading digits to 10 and Trailing digits to 20, you get

> 2^30000;

7940903519[...9001 digits...]35833790694215909376

The problem is apparently one involving simplification: at some point Maple doesn't recognize that sqrt(6)=sqrt(3)*sqrt(2).  If you replace sqrt(6) by sqrt(2)*sqrt(3), you get the correct Jordan form.  Note also:

> Q:= JordanForm(A, output='Q');

  which is the Matrix that should produce Q^(-1) . A . Q = J.  But in fact, this Q is singular: its second column simplifies to all 0.

> simplify(Q);
                   [       1                    2        ]
                   [       -         0          -        ]
                   [       3                    3        ]
                   [                                     ]
                   [    1  (1/2)            1  (1/2)     ]
                   [  - - 2          0      - 2          ]
                   [    6                   6            ]
                   [                                     ]
                   [1  (1/2)  (1/2)       1  (1/2)  (1/2)]
                   [- 2      3       0  - - 2      3     ]
                   [6                     6              ]

You can "extract" a function from a numeric solution something like this (I'll try a rather less trivial example):

> sys:= {diff(x(t),t) = t*y(t) + x(t)*y(t), diff(y(t),t) = y(t) - x(t), x(0) = 1, y(0) = 0};
  Sol:= dsolve(sys, numeric, output=listprocedure);
  X:= subs(Sol, x(t)); Y:= subs(Sol, y(t));

Then X and Y are procedures that will return numerical values: e.g.

> X(Pi); 

−2.89145099805825588

You can also plot them.

However, you can't differentiate them.  What you can do, though, is use the system of differential equations to express their derivatives in terms of X and Y.  Thus:

> Xp:= subs(x=X,y=Y, subs(sys, diff(x(t),t)));

And then this can be plotted:

> plot(Xp, t = 0 .. 5);

These three equations (with m=1) have no real solutions, so it's not surprising that fsolve can't find any.  However, you might try m = 1.9.

>   fsolve(eval(eqs,m=1.9), {x1 = 0 .. 1, x2 = 0 .. 1, x3 = 0 .. 1});

{x1 = .1320826971, x2 = .9721759684, x3 = .7957413345}

This doesn't have x1,x2 and x3 in the right order, but by symmetry you can permute the values so that x3 <= x2 <= x1.

 For the plot, you could try something like this (but it won't give you the interval of m values you mentioned):

> pts:= NULL: 
  for X3 from 0 to 1 by 0.01 do
    R:= fsolve({-3*x1+4*x1^3-3*x2+4*x2^3-3*X3+4*X3^3 = 0, 
         5*x1-20*x1^3+16*x1^5+5*x2-20*x2^3+16*x2^5+5*X3-20*X3^3+16*X3^5 = 0},
        {x1=X3..1,x2=X3..1});
    if type(R,set) then 
      X2:= min(map(rhs,R));
      X1:= max(map(rhs,R));   
      M1:= X1 + X2 + X3;
      pts:= pts, [X1,M1],[X2,M1],[X3,M1];
    end if
  end do:
  plots[pointplot]([pts]);      

For example:

    |\^/|     Maple 13 (IBM INTEL NT)
._|\|   |/|_. Copyright (c) Maplesoft, a division of Waterloo Maple Inc. 2009
 \  MAPLE  /  All rights reserved. Maple is a trademark of
 <____ ____>  Waterloo Maple Inc.
      |       Type ? for help.
> currentdir(); # this returns the current directory
                        "C:\Program Files\Maple 13.02"

> currentdir("d:/test");  # this changes the directory to "d:/test", but returns the old current directory
                        "C:\Program Files\Maple 13.02"

> currentdir(); # check that it worked
                                   "d:\test"

> read "myfile.mpl"; # should read it from d:\test

Note that you use "/" as the directory separator, not "\": you could also use "\\".

Your equations are second order in both v[ds] and v[gs], so you need two more initial conditions to determine the solution uniquely. That's why there are still two arbitrary constants in the solution. 

> eqns := 1.00011000*10^(-7)*(diff(v[gs](t), t))-1.000*10^(-12)
    *(diff(v[ds](t), t))+v[gs](t)+1.000*10^(-18)*(diff(v[gs](t), t, t))
    +5.00*10^(-19)*(diff(v[ds](t), t, t)) = 12, 2.900*10^(-18)*
    (diff(v[ds](t), t, t))+5.00*10^(-7)*(diff(v[gs](t), t))+6.00*10^(-19)
    *(diff(v[gs](t), t, t))+v[ds](t) = 24;
  ics:= v[gs](0) = 3.5, 
      v[ds](0) = `#msub(mi("V"),mo("in",fontweight = "bold"))`,
      D(v[gs])(0)= a, D(v[ds])(0)=b;

  dsolve({eqns, ics});

 To evaluate the RootOfs numerically, you can use evalf.

> evalf(%);

  The result is rather complicated, but that's life.