Just expand the exponentials and use the rule for the cosine:

exp2cos:=A::algebraic*exp(I*B::algebraic)+A::algebraic*1/exp(I*B::algebraic)=2*A*cos(B):
applyrule(exp2cos,expand(r));

                                    2      2
                   1/2   1/2       k  sigma
            sigma 2    Pi    exp(- ---------) cos(k theta)
                                       2

Or, without expanding:

exp2cosh:=A::algebraic*exp(a::algebraic)+A::algebraic*exp(b::algebraic)=
2*A*cosh((a-b)/2)*exp((a+b)/2):

applyrule(exp2cosh,r);

                                    2      2
                   1/2   1/2       k  sigma
            sigma 2    Pi    exp(- ---------) cos(k theta)
                                       2

If the last one is E(x)*E(y) = -1/n if x<>y and E(x)^2 = 1/n, you can apply such rules by applyrule as in these examples:

r:=[E(x::name)*E(y::name) = -1/n,E(x::name)^2 = 1/n]:
applyrule(r,E(u1)*E(u2));
                                1
                              - -
                                n

map2(applyrule,r,[E(u1)*E(u2)*E(u3)*E(u4),E(u1)^2*E(u2)*E(u3),E(u1)^2*E(u3)^2]);
                         [1     1   1 ]
                         [--, - --, --]
                         [ 2     2   2]
                         [n     n   n ]

For expressions similar to u(x)*u(y), the following procedure will do it, by expanding F(u+v)-F(u), and using rules to remove nonlinear terms in v:

fd:=proc(F,u::name,v::name)
local r:
r:=[A::algebraic*v(a::name)*v(b::name)=0,conditional(A::algebraic*v(b::name)^n::posint=0,n>1)]:
subs(u=u+v,F)-F;
simplify(%);
applyrule(r,%);
end proc:

fd(u(x)*u(y),u,v);
                        u(x) v(y) + v(x) u(y)

fd(u(x)*u(y)^2,u,v);
                                                2
                    2 u(x) u(y) v(y) + v(x) u(y)

fd(u(x)*u(y)*u(z),u,v);
           u(x) u(y) v(z) + u(x) v(y) u(z) + v(x) u(y) u(z)

Tracing shows that int(sqrt(1 + sec(x)^2), x) is transformed by the change of variable u=tan(x/2), producing an integrand with a fourth power:

[...]
{--> enter int/algebraic2/algebraic, args = 1/(_X^2-1)*(2*_X^4+2)^(1/2)/(1+_X^2)
[...]

and probably, because of it, it is dispatched to the elliptic integration routines. On the other hand, int(sqrt(2 + tan(x)^2), x) produces a quadratic term:

[...]
{--> enter int/algebraic2/algebraic, args = (2+_X^2)^(1/2)/(1+_X^2)
[...]

and it is integrated in terms of trigonometrics.

First, for such symbolic computation it is advisable to use the inert form of the differential functions like %diff and %Laplacian, instead of the active forms diff and Laplacian that you would use for explicit calculations. Then, the required transformation rules, like permuting the order of the derivatives, can be performed by applyrule:

e2 := %diff(u(x, y), x)-%diff(v(x, y), y);
G2:=map(%Laplacian,e2):
Sxx := map(%diff,G2, x):
Syy := map(%diff,G2, y):
r1:=%Laplacian(%diff(y::function,x::name))=%diff(%Laplacian(y),x):
r2:=%Laplacian(-1*A::algebraic)=-%Laplacian(A):
r3:=%diff(-1*A::algebraic,x::name)=-%diff(A,x):
s:=%Laplacian(u(x,y))=P,%Laplacian(v(x,y))=Q:
applyrule([r1,r2,r3,s],Sxx);

                         / 2   \   /  2    \
                         |d    |   | d     |
                         |--- P| - |----- Q|
                         |  2  |   \dx dy  /
                         \dx   /

applyrule([r1,r2,r3,s],Syy);

                         /  2    \   / 2   \
                         | d     |   |d    |
                         |----- P| - |--- Q|
                         \dy dx  /   |  2  |
                                     \dy   /

Note that linearity of the inert differential functions has to be imposed explicitly, e.g. by using map. And no package is needed for these symbolic computations.

Several steps are involved for arriving to the answer, the same as you would do by hand. So, the main issue is implementing every of these steps as a rule. Also, as these symbolic computations do not make use the inner computational routines associated with function calls like sum(...) or int(...), it is more efficient avoiding them altogether by using the inert forms like Sum(...) or Int(...):

eq := I*hbar*(Sum((diff(c[n](t), t))*f[n](r)*exp(-I*omega[n]*t), n = l .. k)) = 
(1/2)*E[0](e_.r_)*e*(Sum(c[n](t)*f[n](r)*omega[n]*(exp(I*t*
(-omega[n]+Omega))+exp(-I*t*(omega[n]+Omega))), n = l .. k)):
req := Int(conjugate(f[m](r))*rhs(eq), r):
leq := Int(conjugate(f[m](r))*lhs(eq), r):

#rules
sumint:=Int(A::algebraic*Sum(v::algebraic,r::equation),u::name)=Sum(Int(A*v,u),r):
intsum:=Sum(v::algebraic*Int(A::algebraic,u::name),r::equation)=Int(Sum(A*v,r),u):
sup:= Int(A::algebraic*conjugate(f[m](r))*f[n](r), r) = A*delta[m, n]:
kron:=[Sum(F::Not(anyindex(identical(n)))*a::anyindex(identical(n))*delta[m,n],r::equation)
=F*a,Sum(F::Not(anyindex(identical(n)))*delta[m,n],r::equation)=F]:

#lhs
applyrule([sumint,sup],Int(conjugate(f[m](r))*leq,r)):
applyrule(intsum,%):
subs(n=m,applyrule(kron,%)):
LHS:=IntegrationTools:-Expand(%):

#rhs
applyrule([sumint,sup],Int(conjugate(f[m](r))*req,r)):
applyrule(intsum,%):
subs(n=m,applyrule(kron,%)):
IntegrationTools:-Expand(%):
RHS:=collect(%,[Int,E[0],c[m]],simplify):

(LHS=RHS)/op(indets(LHS,specfunc(anything,Int)));

       /d         \
  hbar |-- c[m](t)| I
       \dt        /
  ------------------- = 1/2 e omega[m]
   exp(omega[m] t I)

        (exp(t (-omega[m] + Omega) I) + exp(-I t (omega[m] + Omega)))

        c[m](t) E[0](e_ . r_)

@Jeff 

Actually, an error can be raised for trying to assign 0, by a typed assignment with assertion failure:

restart:
kernelopts(assertlevel=2):

a::Non(0):=0;
Error, assertion failed in assignment, expected Non(0), got 0

a;
                                 a

I would do it in two stages, one for the function call using anyfunc, and then for its arguments, like:

g(u+v*w): 
patmatch(%,f::anyfunc(`+`(algebraic)),'p'):p;
patmatch(op(subs(%,f)),a::name*b::name +c::name,'p'):p;

                           [f = g(u + v w)]
                        [a = v, b = w, c = u]

Note that the type `+`(algebraic) is simple but a bit more general than the pattern in the next line. A more strict (but longer) type version is:

g(u+v*w): 
patmatch(%,f::anyfunc(`&+`(name,`&*`(name,name))),'p'):p;
patmatch(op(subs(%,f)),a::name*b::name +c::name,'p'):p;

                           [f = g(u + v w)]
                        [a = v, b = w, c = u]

Variants may be needed, depending on details. Also, if you do know the name of a variable, v say, and you want linearity in it, you could do instead something like:

g(u+v*w): 
patmatch(%,f::anyfunc(`+`(linear(v))),'p'):p;
patmatch(op(subs(%,f)),a::name*v +c::name,'p'):p;

                           [f = g(u + v w)]
                            [a = w, c = u]

@acer 

I think that such a simple algebraic transformation should be possible to be dealt with also simply, in uniform way, without much need of creativity or fine tuning, as is the case of using applyrule for the above examples:

csq:=conditional(a::algebraic+2*b::algebraic,a=b^2)=(b+1)^2-1:
applyrule(csq,x^2+2*x);
                                    2
                             (x + 1)  - 1

applyrule(csq,x^4+2*x^2);
                              2     2
                            (x  + 1)  - 1

applyrule(csq,x^4/sin(x)^4+2*x^2/sin(x)^2);

                          /   2       \2
                          |  x        |
                          |------- + 1|  - 1
                          |      2    |
                          \sin(x)     /

Just assuming real:

int((r^2+r*cos(theta))*r*(a^2-r^2)^(1/2),[r = 0 .. a,theta = 0 .. 2*Pi]) assuming real;
                               5
                         4/15 a  Pi signum(a)

See this thread, where problems were also reported when a radiation boundary condition is present.

You can easily make the transformations from exp to cosh and sinh with applyrule:

applyrule~(1/2*exp(a::algebraic)+1/2*exp(b::algebraic)=cosh((a-b)/2)*exp((a+b)/2),p):
applyrule~(1/2*exp(a::algebraic)-1/2*exp(b::algebraic)=sinh((a-b)/2)*exp((a+b)/2),%);

I would use:

ee:=itail*(exp(vd/Vth)-1)/(1+exp(vd/Vth)):
applyrule((exp(z::algebraic)-1)/(exp(z::algebraic)+1)=tanh(z/2),ee);

                                vd
                          tanh(-----) itail
                               2 Vth