You don't have what's usually called an "analytic" solution.  It's a numeric solution. 

Now lnH2 is a function (so that e.g. lnH2(3.5)  should give you a numeric result).  Therefore you might normally use

  evalf(Int(exp(lnH2(x)),x=0..1));

except that lnH2(x) with x as a symbolic variable doesn't work. You could try

 sol := dsolve({initial, syst}, func, numeric, output=listprocedure);
 lhH2:= subs(sol,lnH(x));
evalf(Int(exp(lnH2(x)), x=0..1));

or with your lnH2 you could try

 evalf(Int(exp @ lnH2, 0 .. 1));

You could get it by making this derivative a new dependent variable v in a system of PDE's.

> SYS:= {v(x,t)=rhs(PDE), diff(u(x,t),t) = v(x,t)};
   pds:= pdsolve(SYS,IBC,numeric, time = t, range = 0 .. 5, 
      timestep = 1/2000, spacestep = 1/2000);
   R:=  pds:-value(output = listprocedure);
   ut:= subs(R, v(x,t));
   evalf(Int(ut(x,2),x=0..5,method=_Dexp));

-.2599816919

You might try extracting the f's using indets.  Thus:

> seq(frontend(coeff, [Expr, op(indets(Expr,f(x[i],anything)))]),i=0..2);

The behaviour you describe is bizarre.  I've never seen anything like it occur in any Maple version.  Can you upload a worksheet where this occurs for you, and maybe a screen shot?  The only thing I can think of is that somehow a font has become corrupted.  I don't know how that could even happen.

There is some hope. 

> RatSys:= convert(RedSys, rational);
    # to avoid possible numerical problems with floats
   RS1:= RatSys assuming t>0,t<20,X97(t) - 1/2500*X63(t)>3/200;
    # one of the cases; casesplit can't handle piecewise
   RS1S:= PDEtools[casesplit](RS1);
   dsolve({op(RS1S), X74(0)=0,X97(0)=0,X21(0)=0,D(X21)(0)=1,
     seq((D@@j)(X21)(0)=0,j=2..5)}, numeric);

The main differential-equation solving command in Maple is dsolve. 

> k2:= a*k1/b;
  des:= { diff(y1(t),t) =  (a-y1(t))/T1 - k1*y1(t) + k2*y2(t),
          diff(y2(t),t) =  (b-y2(t))/T2 - k2*y2(t) + k1*y1(t) };
  dsolve(des);

You didn't say what the "(->)" was.  If you use Simplify. Symbolic (I assume it's from the context menu), you get the same thing on both sides.  If you use Simplify, Simplify you don't, because (without any assumptions), the two are not necessarily equivalent.  Maple doesn't assume anything is real unless you tell it to, and identities such as (x*y)^p = x^p * y^p are not always true for complex numbers.  For example, the two sides of your equation are not the same if beta = I and m = -I (where I = sqrt(-1)).

My first idea was to use display to combine several contourplot3d's.  Something like
 

> with(plots):
  expr := x^2 + y^2; V := [1,2,3,4]: C:= [red,yellow,green,blue]:
  display([seq(contourplot3d(expr, x=-2..2, y=-2..2, contours = [V[j]], colour=C[j]),
     j=1..4)]);   

where V is a list of contour values and C a list of colours.  However, that doesn't work: this seems to be a bug.  It would have worked in contourplot.  So you could try something like this:

> with(plots):
  expr := x^2 + y^2;  V := [1,2,3,4]: C:= [red,yellow,green,blue]:
  T := plottools[transform](unapply([x,y,expr],x,y)):
  T(display([seq(contourplot(expr,x=-2..2, y=-2..2, contours=[V[j]], colour=C[j]),
    j=1..4)]));  

listplot in the plots package will plot a Vector (in the form of points [i, V[i]], i = 1 to the dimension of the Vector, joined by straight lines).  If you don't want the points joined, you can use the option style=point.
Several of them can be combined using display.  Thus:

> with(plots):
  display([listplot(X), seq(listplot(U[1..10,j], colour=HUE(j/20), style=point),j=1..20)]);

One way to incorporate a constraint in dfieldplot is to make the direction field non-real when the constraint is violated.  Thus:

> dfieldplot([diff(x(N),N)=-3*x(N)+sqrt(6)*(y(N))^2+3/2*x(N)*((x(N))^2+1-(y(N))^2), 
  diff(y(N),N)=-sqrt(6)*x(N)*y(N)+3/2*y(N)*((x(N))^2+1-(y(N))^2)+Heaviside(x(N)^2+y(N)^2-1)*I],
 [x(N),y(N)],N=0..20, x=-1..1, y=0..1, color=black, dirfield=[30,30]); 

For your second plot, I'd use dsolve and odeplot.  Thus:

> S:= dsolve([diff(x(N),N)=-3*x(N)+sqrt(6)*(y(N))^2+3/2*x(N)*((x(N))^2+1-(y(N))^2), 
       diff(y(N),N)=-sqrt(6)*x(N)*y(N)+3/2*y(N)*((x(N))^2+1-(y(N))^2),
       x(0) = 0, y(0) = 0.01], numeric);
  plots[odeplot](S, [[N,x(N)],[N,x(N)^2 + y(N)^2]], N = 0 .. 10, colour=[red,blue], 
       legend = [x, x^2+y^2],labels=[N,""]);

> add(add(list1[i,j]*list2[i]*list2[j], j=i..3), i=1..3);

Note that the limits in the inner sum can depend on the index variable for the outer sum, but not vice versa.

Your procedure VWN5 uses lexical scoping, i.e. there are many local variables of VWN5 which are used (but not declared) in the procedure eps, which is defined within VWN5.  Actually I don't think codegen[fortran] (or the more recent CodeGeneration[Fortran]) can handle procedures defined within other procedures.
There are other problems with your procedure as well.  For example, b and c are used in Q and the procedure X, but apparently have not been assigned values.
Translation to Fortran (and other languages) only works for procedures that are purely numerical: they can't handle expressions involving unassigned variables.

What have you tried so far?  What point are you stuck on?

Whether one expression is considered as "simpler" than another may be a matter of context or personal taste.  It's not at all clear that tan(x) is simpler than sin(x)/cos(x).  Mathematica may make one choice, Maple another.  But there is certainly room for improvement in the simplify command's power and flexibility.

MathCad?  I hope you mean Maple.

One thing I see right away is that you have

assign(delta, phi)

where you probably want to assign the results of the previous solve command: that would be

assign(%[1])

Another problem: you have  L[pile[i]] where you probably want L[pile][i] (and similarly for other variables)