Are you saying f(e) is quadratic in w?  So let f(e) =   A*w^2 + B*w + C.

> e := ((w-b)/b)^lambda;
  prof:= f(e) - w*L;
  prof1:= eval(prof, f(e) = A*w^2 + B*w + C);

Now you say the derivative of this wrt w is 0 at w = -b/(lambda - 1).

> eval(diff(prof1, w), w = -b/(lambda-1));

That being 0 gives you a relation between A and B, which you can solve for either of these, let's say for B.

> B0:= solve(%, B);

B0 := (2*A*b+L*lambda-L)/(lambda-1)

The maple tag is acting up again: that should be

    B0 := (2*A*b+L*lambda-L)/(lambda-1)

So now f(e) = A*w^2 + B0 * w + C.  Now.we can express w in terms of e.  But since I defined e as an expression in w, I'll use E instead.

> wE:= solve(e=E, w);

    wE := exp(ln(E)/lambda)*b+b

And then the function f is

> f(E) =  simplify(A*wE^2 + B0*wE + C);

f(E) = (A*b^2*E^(2/lambda)*lambda-A*b^2*E^(2/lambda)+2*A*b^2*E^(1/lambda)*lambda+A*b^2*lambda+A*b^2+b*L*lambda*E^(1/lambda)+b*L*lambda-b*L*E^(1/lambda)-b*L+C*lambda-C)/(lambda-1)

Is that what you had in mind?

In general, Maple ignores the possibility of "special" values of variables.  Otherwise it could never divide by anything without being sure that the denominator is nonzero.  So for example int(x^n, x) = x^(n+1)/(n+1), although that is invalid for the special value n = -1.  This usually works out OK, but if you want to evaluate a general expression at a particular value of a variable you have to worry that this value might be one of those "special" values for which the expression is invalid. 

In general, Maple ignores the possibility of "special" values of variables.  Otherwise it could never divide by anything without being sure that the denominator is nonzero.  So for example int(x^n, x) = x^(n+1)/(n+1), although that is invalid for the special value n = -1.  This usually works out OK, but if you want to evaluate a general expression at a particular value of a variable you have to worry that this value might be one of those "special" values for which the expression is invalid. 

If w = exp(I*theta) = cos(theta) + I*sin(theta), exp(I*theta/3) = cos(theta/3) + I*sin(theta/3) is one of the three cube roots of exp(I*theta).  That is, if
alpha = exp(2*Pi*I/3), exp(I*theta/3) is either w^(1/3), alpha*w^(1/3) or alpha^2*w^(1/3), and then cos(theta/3) = (exp(I*theta/3) + 1/exp(I*theta/3))/2 is one of

(w^(1/3) + 1/w^(1/3))/2;

(alpha*w^(1/3) + 1/(alpha*w^(1/3)))/2

(alpha^2*w^(1/3) + 1/(alpha^2*w^(1/3)))/2

Which one it is will vary: for example, it is the first for -Pi <= theta <= Pi, the second for
Pi < theta <= 3*Pi, the third for 3*Pi < theta <= 5*Pi.

Note also that there is no expression for cos(theta/3) in terms of "real radicals":  you must use cube roots of complex quantities.

> w:= exp(I*theta);
  alpha:= exp(2*Pi*I/3);
  v:= w^(1/3);
  r1:= (v+1/v)/2; r2:= (alpha*v + 1/(alpha*v))/2;
  r3:= (alpha^2*v + 1/(alpha^2*v))/2;
  plot([r1, r2, r3], theta = -Pi .. 5*Pi, 
     colour=[red,blue,green]);

If w = exp(I*theta) = cos(theta) + I*sin(theta), exp(I*theta/3) = cos(theta/3) + I*sin(theta/3) is one of the three cube roots of exp(I*theta).  That is, if
alpha = exp(2*Pi*I/3), exp(I*theta/3) is either w^(1/3), alpha*w^(1/3) or alpha^2*w^(1/3), and then cos(theta/3) = (exp(I*theta/3) + 1/exp(I*theta/3))/2 is one of

(w^(1/3) + 1/w^(1/3))/2;

(alpha*w^(1/3) + 1/(alpha*w^(1/3)))/2

(alpha^2*w^(1/3) + 1/(alpha^2*w^(1/3)))/2

Which one it is will vary: for example, it is the first for -Pi <= theta <= Pi, the second for
Pi < theta <= 3*Pi, the third for 3*Pi < theta <= 5*Pi.

Note also that there is no expression for cos(theta/3) in terms of "real radicals":  you must use cube roots of complex quantities.

> w:= exp(I*theta);
  alpha:= exp(2*Pi*I/3);
  v:= w^(1/3);
  r1:= (v+1/v)/2; r2:= (alpha*v + 1/(alpha*v))/2;
  r3:= (alpha^2*v + 1/(alpha^2*v))/2;
  plot([r1, r2, r3], theta = -Pi .. 5*Pi, 
     colour=[red,blue,green]);

Try S[1,1].

Notice that there is a difference between a 5-element column Vector and a 5 x 1 Matrix.  S[1] would work if S was a Vector.

Try S[1,1].

Notice that there is a difference between a 5-element column Vector and a 5 x 1 Matrix.  S[1] would work if S was a Vector.

I suspect that what JMoin wants to plot is not the triple integral itself but the region of three-dimensional space being integrated over.  Thus for an integral such as

> Int(Int(Int(f(x,y,z), x = -z .. y), y = 0 .. z), z = 0 .. 1);

the question is what does the region look like?

Maple doesn't plot 3-dimensional regions directly, but it does plot surfaces in three dimensions. 
One way is to plot each face of the region separately, something like this:

> with(plots):
  P1:= plot3d([-z,y,z],y=0..z,z=0..1,colour=red):
  P2:= plot3d([y,y,z],y=0..z,z=0..1,colour=red):
  P3:= plot3d([x,0,z],x=-z..0,z=0..1,colour=green):
  P4:= plot3d([x,z,z],x=-z..z,z=0..1,colour=green):
  P5:= plot3d([x,y,1],x=-1..y,y=0..1,colour=blue):
  display(P1,P2,P3,P4,P5,axes=box,scaling=constrained,
     style=patchnogrid,labels=[x,y,z], orientation=[35,75]);

Another way is to consider the inequalities -z <= x, x <= y, 0 <= y, y <= z, 0 <= z and z <= 1,
and write each in the form something <= 0.  You want all the left sides <=0, which is equivalent to
saying the maximum of the left sides <= 0.  So:

> F := max (-z - x, x - y, -y, y - z, -z, z - 1);   
  plots[implicitplot3d](F, x = -1.1 .. 1.1, y = -0.1 .. 1.1, z = -0.1 .. 1.1,
       axes = box, labels = [x,y,z], scaling = constrained, style=patchnogrid, 
       orientation = [35,75], lightmodel = light1); 

I suspect that what JMoin wants to plot is not the triple integral itself but the region of three-dimensional space being integrated over.  Thus for an integral such as

> Int(Int(Int(f(x,y,z), x = -z .. y), y = 0 .. z), z = 0 .. 1);

the question is what does the region look like?

Maple doesn't plot 3-dimensional regions directly, but it does plot surfaces in three dimensions. 
One way is to plot each face of the region separately, something like this:

> with(plots):
  P1:= plot3d([-z,y,z],y=0..z,z=0..1,colour=red):
  P2:= plot3d([y,y,z],y=0..z,z=0..1,colour=red):
  P3:= plot3d([x,0,z],x=-z..0,z=0..1,colour=green):
  P4:= plot3d([x,z,z],x=-z..z,z=0..1,colour=green):
  P5:= plot3d([x,y,1],x=-1..y,y=0..1,colour=blue):
  display(P1,P2,P3,P4,P5,axes=box,scaling=constrained,
     style=patchnogrid,labels=[x,y,z], orientation=[35,75]);

Another way is to consider the inequalities -z <= x, x <= y, 0 <= y, y <= z, 0 <= z and z <= 1,
and write each in the form something <= 0.  You want all the left sides <=0, which is equivalent to
saying the maximum of the left sides <= 0.  So:

> F := max (-z - x, x - y, -y, y - z, -z, z - 1);   
  plots[implicitplot3d](F, x = -1.1 .. 1.1, y = -0.1 .. 1.1, z = -0.1 .. 1.1,
       axes = box, labels = [x,y,z], scaling = constrained, style=patchnogrid, 
       orientation = [35,75], lightmodel = light1); 

Once you have the Plot component, you can use DocumentTools[SetProperty]  if you wish to change the
size under programmatic control, rather than by clicking.  For example, if the Plot component is Plot0
(which will be the name of the first Plot component you insert in your worksheet):

> use DocumentTools in
    SetProperty(Plot0, 'PixelHeight',250);
    SetProperty(Plot0, 'PixelWidth', 250);
    SetProperty(Plot0, 'value', plot3d(x^2-y^2,x=-1..1,y=-1..1))
end use:

Note that components are only available in the Standard GUI, not in Classic.

Once you have the Plot component, you can use DocumentTools[SetProperty]  if you wish to change the
size under programmatic control, rather than by clicking.  For example, if the Plot component is Plot0
(which will be the name of the first Plot component you insert in your worksheet):

> use DocumentTools in
    SetProperty(Plot0, 'PixelHeight',250);
    SetProperty(Plot0, 'PixelWidth', 250);
    SetProperty(Plot0, 'value', plot3d(x^2-y^2,x=-1..1,y=-1..1))
end use:

Note that components are only available in the Standard GUI, not in Classic.

Actually you didn't post the corrected code.  Your code had e^... instead of exp, and a lot of what appear to be function calls but may be intended to be multiplications.  Are things such as `&lambda;bj` intended to be actual names of parameters or do you mean lambda*b*j, or what?  You have some sums over i and j, but also some j's appearing outside of any sums.  

If you want symbolic solutions, you should use solve.  However, it's unlikely to produce any results on complicated non-polynomial systems.  Or if you want a numerical solution (when all parameters are given numerical values) you can try fsolve.

Actually you didn't post the corrected code.  Your code had e^... instead of exp, and a lot of what appear to be function calls but may be intended to be multiplications.  Are things such as `&lambda;bj` intended to be actual names of parameters or do you mean lambda*b*j, or what?  You have some sums over i and j, but also some j's appearing outside of any sums.  

If you want symbolic solutions, you should use solve.  However, it's unlikely to produce any results on complicated non-polynomial systems.  Or if you want a numerical solution (when all parameters are given numerical values) you can try fsolve.

I assume your system is linear, otherwise such an expansion will not generally exist.  You can get a fundamental set of solutions of the homogeneous problem by using different initial conditions (one variable = 1, the others = 0).

For example:

> sys:= {D(x)(t) = 2*x(t) - t*y(t) + cos(t), 
         D(y)(t) = x(t) - t*z(t), 
         D(z)(t) = y(t) - t*z(t)};
   homog:=  {D(x)(t) = 2*x(t) - t*y(t) , 
             D(y)(t) = x(t) - t*z(t), 
             D(z)(t) = y(t) - t*z(t)};
  S1:= dsolve(homog union {x(0)=1, y(0)=0, z(0)=0}, numeric, output=listprocedure);
  S2:= dsolve(homog union {x(0)=0, y(0)=1, z(0)=0}, numeric, output = listprocedure);
  S3:= dsolve(homog union {x(0)=0, y(0)=0, z(0)=1}, numeric, output = listprocedure);
  S4:= dsolve(sys union {x(0)=0, y(0)=0, z(0)=0}, numeric, output = listprocedure);
 

Now the general solution will be  x(0)*S1 + y(0)*S2 + z(0)*S3 + S4. Thus the coefficients in x(t)
would be given by
 

> e1:= subs(S1, x(t));
  e2:= subs(S2, x(t));
  e3:= subs(S3, x(t));
  e4:= subs(S4, x(t));

Start with the following command:

> plots[setoptions](axesfont = [Courier, 20]);

... and then call up the Curve Fitting Assistant; it should use 20-point Courier for the numbers on the axes.