I hope I understand what you're trying to do here: you want to give makeplot the name of a variable to which a plot is assigned, and you want makeplot to produce the associated eps file using that name with either _classic or _standard inserted.  Something like this?

> makeplot:= proc(p::evaln)
     local fn, opts;
     if IsWorksheetInterface('Standard') then
        fn:= cat(convert(p,string), "_standard.eps");
        opts:= "color,portrait,noborder"
     else
         fn:= cat(convert(p,string), "_classic.eps");
         opts:= "nocolor,portrait,noborder"
     end if;
     plotsetup(ps, plotoutput=fn, plotoptions=opts);
     print(eval(p));
     plotsetup(default);
   end proc;

Then you can try e.g.

> for i from 1 to 5 do
        P[i] := plot(cos(i*x), x=-2*Pi..2*Pi);
        makeplot(P[i])
    end do:

solve gives a solution involving a RootOf.  To get both solutions, use allvalues.

> solve({E1,E2},{x,y});

{x = 30*RootOf(11*_Z^2+8*_Z+1,label = _L4), y = 60*RootOf(11*_Z^2+8*_Z+1,label = _L4)+40}

> allvalues(%);

{x = -120/11+30/11*5^(1/2), y = 200/11+60/11*5^(1/2)}, {x = -120/11-30/11*5^(1/2), y = 200/11-60/11*5^(1/2)}

> evalf(%);

{x = -4.810723701, y = 30.37855260}, {x = -17.00745812, y = 5.98508376}

You seem to be stuck on variations of this one problem.  From what I told you in another thread, the only way to have y(infinity) = 0 in a solution of diff(y(r), r, r)+2*(diff(y(r), r))/r = f(y) with f continuous at 0 would be f(0) = 0.  Now here you have f(y) = - C*polylog(3/2, -exp(C2*(C3-y))).  But polylog(3/2, t) has no real zeros other than 0, so except in the trivial case C=0 your boundary condition is impossible.

You can actually get symbolic solutions of a sort:

> des:= {diff(B1(z),z)= - a*B1(z)-(g/AE)* B4(z)* B1(z),
 diff(B2(z),z)= +a*B2(z) + (g/AE)* B3(z)* B2(z),
 diff(B3(z),z)= - a*B3(z) + (g/AE)* B2(z)* B3(z),
 diff(B4(z),z)= + a*B4(z)-(g/AE)* B1(z)* B4(z)};
> dsolve(des);

But I don't know how much that will help you.
One thing that might be useful is to non-dimensionalize your equations: by scaling the dependent and independent variables, you can get rid of two parameters, e.g. reducing to the case a=1, g/AE=1.

For numerical values of the parameters, you might use dsolve(..., numeric).
Thus:

> des1:= eval(des, {a=1,g=AE});
   bcs1:= {B1(0)=3,B2(2)=4*B1(2), B3(2)=3,B4(2)=2*B3(2)};
   Sol:= dsolve(des1 union bcs1, numeric);
   plots[odeplot](Sol, [[z,B1(z)],[z,B2(z)],[z,B3(z)],[z,B4(z)]],
      legend=[B1,B2,B3,B4]);

The problem is that absolute value is a non-differentiable function, and the Optimization procedures don't work well with that.  But in many cases you can avoid the use of the absolute value by using constraints.  For example, to minimize |X| (where X is some expression), you could minimize t subject to constraints t >= X and t >= -X.

> max(indets(MM,integer));

ImageTools:-Preview produces a 3d plot, and you probably want a 2d plot. It is possible to make an image into a 2d plot.  It's simplest with a black-and-white image.

> FileLocation:= cat(kernelopts(datadir), "/images/phone.jpg");
   A:= ImageTools:-Read(FileLocation):
   P:= plots[listdensityplot](ImageTools:-Rotate(A,90), 
      style=patchnogrid):

The result can be used as background for an animation.  Warning: this tends to produce huge plot objects that the worksheet interface may find hard to handle if the number of frames is not very small.

> plots[animate](plot, [t*x, x=0..255], t=0..1, background=P, 
    frames=4);

In this case the error message is quite correct.  Calling A(5,5) will produce a huge recursion, with very many levels (i.e. function calls active at the same time).  The Ackermann function is an important example in the theory of recursive functions, but it's not something that's practical to actually compute, except for very small inputs.

Something like this, perhaps (producing a list rather than an array):

> H := [seq]( `if`(pointsd[2][j+1]<pointsd[2][j], j-1, NULL), j=1..T);

Another way:

> H := map(t -> t-1, select(t -> (pointsd[2][j+1]<pointsd[2][j]), [$1..T]));

There may not be closed-form solutions, which is what solve looks for, except for the trivial: x = (+/-) sqrt(a) for v > 0.  For numerical values of the parameters, you can use fsolve to find a single real numerical solution, or RootFinding[Analytic] for complex solutions.  For example:

> RootFinding[Analytic](BesselK(1, x)*x-(1-x^2)^(1/2)*BesselK(0, x)*BesselJ(1, (1-x^2)^(1/2)), x = -5 - 5*I .. 5 + 5*I);

1.763756716-3.383224668*I, .1190027667, 1.763756716+3.383224668*I

I think you want something like this:

> for i from 1740 to 1770 do
     Y[i]:= readdata(cat("/Users/humadih/test/",i,".tmp"), 1)
   end do;

This should produce separate lists Y[1740] to Y[1770].

If you want a Matrix with these as rows, you can then do

> Matrix([seq(Y[i],i=1740..1770)]);

The reason is that (f[x]-z)/(f(x)-y) doesn't occur as a syntactic unit in A.  There's a complication, apparently caused by the use of x as a subscript as well as a variable, but this works:

> A := y*(f[X]-z)/(f(x)-y);
  simplify(simplify(A,{(f[X]-z)/(f(x)-y)=lambda}));

lambda*y

s and k are nonnegative integers, I assume.

> J:= Int(x^s * exp(-x*(1+k)/2 - exp(-x)/2), x=-infinity .. infinity);
  F:=  IntegrationTools[Change](J,x=-ln(t));

 F := Int((-ln(t))^s*t^(-1/2+1/2*k)*exp(-1/2*t),t = 0 .. infinity)

Now the case s = 0 can be done in closed form:

> F0:= value(eval(F,s=0));

F0 := 2^(1/2+1/2*k)*GAMMA(1/2+1/2*k)

And the general case is (-2)^s times the s'th derivative of that with respect to k.

> (-2)^s * diff(F0, k$s);

(-2)^s*Sum(binomial(s,_k1)*diff(2^(1/2+1/2*k),`$`(k,_k1))*diff(GAMMA(1/2+1/2*k),`$`(k,s-_k1)),_k1 = 0 .. s)

Usually this does indicate a singularity, and indeed this one does have singularities. 
What initial condition are you using?

The solution of this equation with initial condition u(t0) = u0 can be written as

t - t0 = int(1/(s^4 - s + 1), s = u0 .. u)

Note that s^4 - s + 1 > 0 for all real s, and int(1/(s^4-s+1), s=-infinity .. infinity) < infinity.    So u -> infinity at

t = t0 + int(1/(s^4-s+1), s=u0 .. infinity).   For example, for initial condition u(0) = 0:

> evalf(Int(1/(s^4-s+1), s = 0 .. infinity);

      1.847670652

> S:= dsolve({diff(u(t),t) = u(t)^4 - u(t) + 1, u(0)=0}, numeric);

      S := proc(x_rkf45) ... end proc

> S(1.9);

Error, (in S) cannot evaluate the solution further right of 1.8476704, probably a singularity