I think this will be problematic due to the presence of the factor exp(-(sin(phi)^2)) in the integrand.

You should be able to integrate it numerically for given values of x, sigma, theta, k and Q

I might be wrong but I have the feeling that what you have done here is not what you intended. What are you trying to achieve ? 

Regarding subs, the error occurs because you are asking Maple to sustitute x<3 into f, which doesn't make sense. It needs to be something=somethingelse. So, subs(p = 1/3, x = 3, f) is OK. Check the help page ?subs

The underlying problem seems to be that GG is unassigned when this line is executed.
GG := [op(GG), ......

So, if you assign the null list to GG when you do the same for G that should solve the problem:
G := []: GG := []:

Edit: The reason, I think, Maple 14 gives "too many levels of recursion" is that when GG is unassigned, op(GG) just returns, GG, so GG:=[op(GG),.... becomes an infinite loop. Ths seems to be a more helpful error on Maple 14 than in Maple 15 where is just loses the kernel connection. I'm using Maple 15 and it took me a while to find the problem due to this.

It may help to use Int instead of int

if alpha = x-2 then why do you substitute x = alpha-2 ? Likewise for the other one. If those are typos then I suppose it highlights the problem with doing the substitutions in that way - you have manually manipulated those expressions for the eval and subs. One way to avoid that is by

restart:
ODE := diff(y(x), x) = (y(x)+1)/(x+2)-exp((y(x)+1)/(x+2)):
par1:=alpha = x+2:
par2:=y(alpha) =  y(x)+1;

par1a:=x=solve(par1,x);par2a:=y(x)=solve(par2,y(x));

                         x = -2 + alpha
                      y(x) = -1 + y(alpha)

par2d:=diff(lhs(par2a),x)=diff(rhs(par2a),alpha);

                   d            d            
                  --- y(x) = ------- y(alpha)
                   dx         dalpha         

subODE:=subs(par2d,lhs(ODE))=eval(rhs(ODE),{par1a,par2a});

             d               y(alpha)      /y(alpha)\
          ------- y(alpha) = -------- - exp|--------|
           dalpha             alpha        \ alpha  /

dsolve(subODE);
             y(alpha) = -ln(ln(alpha) + _C1) alpha

Or I maybe I misunderstood what you are doing ?

Supposing you import the data as an array, and the array is numeric, this plots the 2nd column (B) vs the third column (C). If the array is called N

plot(N[..,2],N[..,3]);

If the first row contains a header (eg "Temp" etc) then you need to remove that. One way to do that is to convert the array to a matrix  with convert(M,Matrix) and then LinearAlgebra:-DeleteRow then you can use the same command to plot.

To add columns, just make a new plot with the new column and use plots:-display to display them. eg

p1:=plot(N[..,2],N[..,3]):
p2:=plot(N[..,1],N[..,3]):
plots:-display(p1,p2);

Your system has no unique solution.

You have 3 equations in 4 unknowns. If you have 4 equations in 4 unknowns (or 3 equations in 3 unknowns), you may get a result your desire if the system is consisent. For example:

sys:={ a + b = 1, c = a + b +1, d = a + c, b = d + 1}:
solve(sys);

                 {a = -1, b = 2, c = 2, d = 1}

Edit: Of course, if you use the 2nd system in your post, you get the result:

sys:={ a +b = 1, c = a +b +1 }:
solve(sys);
                   {a = 1 - b, b = b, c = 2}

Have you tried resizing the plot in Maple ?

This is what I see on my machine, in the attached screencap The arrows look OK to me.

.

As for exporting, I haven't used maple for a while for presentation, but the last time I did, I used rtf format, then edited it in MS Word

Edit: Oh, you mean the export format for the graphic not the worksheet. I don't have much experience with that, but I always found bmp was easy too work with and should be lossless.

The problem seems to be caused by this line

PVBC := int((p-p*delta)(SalesBCM)/(1+r)^t, t = 0 .. cmtime)

I believe it should be 

PVBC := int((p-p*delta)*(SalesBCM)/(1+r)^t, t = 0 .. cmtime)

In your worksheet, this evaluated to

PVBC := int((50-50*delta((50+50*delta)*(0.3e-1+(.36*(1-exp(-.39*t)))/(12.00000000*exp(-.39*t)+1))*(1-(1-exp(-.39*t))/(12.00000000*exp(-.39*t)+1))))/1.08^t, t = 0 .. -2.564102564*ln((1.+5.*delta)/(53.+5.*delta)))

This resulted in delta((....  being treated as the function delta(.....) rather than delta*((......

This doesn't become apparent until you try to find the stationary points by differentiation, whereupon it becomes clear that the derivatives are not as expected. 

Edit:typo

To extract a Row from a matrix you can use Row For example:

Download Row.mw

Without seeing the rest of your code, I'm not sure, but have you tried using

evalf(fsolve(xi_M[Temp,sig,interfaceenergy)=0.5,T, T=100..400]))

?

It may be preferable if you could post your code, (after removing the output of the worksheet if long expressions are involved)

Regarding the 2nd question, if F has a normal distribution truncated at [a,b] , it's PDF is 

pdf_tn:=PDF(TruncNormal(mu,sigma,a,b),x);

         /                   /          /          2\\//                
         |                   | (1/2)    |  (x - mu) || |    (1/2)       
piecewise|a          |                   |          |         2 || \                
         \                   \          \  2 sigma  //                  

  /     /          (1/2)\        /          (1/2)\\\   \
  |1    |(b - mu) 2     |   1    |(a - mu) 2     |||   |
  |- erf|---------------| - - erf|---------------|||, 0|
  \2    \    2 sigma    /   2    \    2 sigma    ///   |
                                                       /

Then, the PDF of a normal distribution with mean mu+t, truncated at [a,b] is

subs(x=y-t,pdf_tn);
         /                           /          /              2\
         |                           | (1/2)    |  (y - t - mu) |
piecewise|a          |                           |          |           2   |
         \                           \          \    2 sigma    /

  \//                /     /          (1/2)\
  | |    (1/2)       |1    |(b - mu) 2     |
  | |2 Pi      sigma |- erf|---------------|
  | \                \2    \    2 sigma    /
  /                                         

          /          (1/2)\\\   \
     1    |(a - mu) 2     |||   |
   - - erf|---------------|||, 0|
     2    \    2 sigma    ///   |
                                /

Note that, formally, to compute this PDF, we need to divide by abs(diff(y(x),x)) 
which is just 1 in this case, so I omitted it.

There was a very interesting discussion on this topic, which included the method shown by epostma above,  recently here:

http://www.mapleprimes.com/questions/123476-Seeking-A-Good-Praxis-In-Solving-Equations

A good starting point could be here:

http://en.wikipedia.org/wiki/Lie_point_symmetry

Since the x range is small, you could use a bivariate taylor expansion of Q. Of course, if you are interested in t>>0  this won't be much use. I got a result, expanding to 8th order around x=0.015 and t=10 and used that in your integral to get a result. I think you would need to look at the taylor polynomial to see how well it resembles Q over the period you are interested in and change the order and t location of the expansion as appropriate.

Download EhsanKayo_02.mw