You might try this.

> pprintf:= proc(S::string)
  local X,L,R,n,i,j,k;
  uses StringTools, Typesetting;
  L:= [RegSplit("%a",S)];
  X:= [args[2..-1]];
  R:= NULL;
  j:= 1;
  for i from 1 to nops(L) do
    n:= CountCharacterOccurrences(L[i],"%");
    R:= R, mi(sprintf(L[i],seq(X[k],k=j..j+n-1)));
    j:= j+n;
    if j <= nops(X) then 
       R:= R, op(map(Typeset,[X[j]]));
       j:= j+1;
    end if;
  end do;
  mrow(R)
end proc;

For example:

> pprintf("Here are integer %d and float %3.2f\n and expression %a", 32, evalf(Pi),Pi/(1+x^2) );

It's hard to answer your question without more details.  Do the examples in the help page ?DynamicSystems,BodePlot work for you?  Could you post the actual commands you're using?

I don't know if it's possible to see what substitution is being made within the IntegrationTutor, but in the worksheet itself you can do something like this:

> with(Student[Calculus1]):
  J := Int(sin(exp(x)) * exp(x), x); 

J := Int(sin(exp(x))*exp(x), x)

> R := Hint(J);

R := [change, u = exp(x), u]

> J := Rule[R](J);

J := Int(sin(exp(x))*exp(x), x) = Int(sin(u), u)

Moreover, in the worksheet there's no problem with infinity as an endpoint.

It seems to me that infnorm(f, x=a..b) initially steps through the interval a..b with steps of  size (b-a)/1000, sometimes taking smaller steps near what appears to be a local maximum.  Moreover, the 1000 is hard-coded: there seems to be no way of adjusting it.

An amusing consequence of this:

> numapprox[infnorm](sin(1000*Pi*x),x=0..1);

0.

test2:= proc(L::uneval)
   assign(map(t -> (t=eval(t)+1), eval(L,1)));
 end proc;

In principle, you can use limit.  For example, f(x) is continuous at x=a if and only if limit(f(x),x=a) = f(a).

You really have just one differential equation (the first one, for U(t)), and two algebraic equations for I1(t) and I2(t) in terms of diff(U(t),t).

> dsolve({sys[1],U(0)=0});

U(t) = Itot*R-exp(-1/R/C*t)*Itot*R

 > eval({sys[2],sys[3]},%);

{1/C*exp(-1/R/C*t)*Itot = (Itot-I1(t))/C, 1/C*exp(-1/R/C*t)*Itot = I2(t)/C}

> solve(%,{I1(t),I2(t)});

{I1(t) = -exp(-1/R/C*t)*Itot+Itot, I2(t) = exp(-1/R/C*t)*Itot}

> f := z^a * sum(z^k/(k+a), k=0..infinity);

f := z^a*LerchPhi(z,1,a)

> R:= solve(f=w,z);

R := RootOf(-_Z^a*LerchPhi(_Z,1,a)+w)

>  S := series(R, w, 4);

This result has the order-0 term RootOf(LerchPhi(_Z,1,a)).   I would guess, however, that (if a > 0) you would want the order-0 term to be 0,

and to get a series in powers of w^(1/a).  Let's call that s, and write the equation as

> eq := z * LerchPhi(z,1,a)^(1/a) = s;

eq := z*LerchPhi(z,1,a)^(1/a) = s

>  series(solve(eq,z), s, 8);

a^(1/a)*s+(-a^(1/a)/(a^(-1/a))/(a+1))*s^2+(-1/2*a^(2/a)*(a^2-a-4)/(a^(-1/a))/(a+2)/(a+1)^2)*s^3+(-1/3*a^(3/a)*(a^4-a^3-10*a^2+4*a+18)/(a^(-1/a))/(a+1)^3/(a+2)/(a+3))*s^4+(-1/24*a^(4/a)*(6*a^7+11*a^6-119*a^5-235*a^4+529*a^3+1016*a^2-632*a-1152)/(a^(-1/a))/(a+1)^4/(a+3)/(a+2)^2/(a+4))*s^5+(-1/10*a^(5/a)*(2*a^9+7*a^8-57*a^7-203*a^6+379*a^5+1378*a^4-888*a^3-3222*a^2+684*a+2400)/(a^(-1/a))/(a+4)/(a+3)/(a+5)/(a+1)^5/(a+2)^2)*s^6+(-1/720*a^(6/a)*(-45078*a^10+9049664*a^3-4942098*a^4-140531*a^9-4936896*a+120*a^13+388306*a^6+252156*a^8+1555497*a^7+159*a^11-6220800+1322*a^12+10065072*a^2-5856333*a^5)/(a^(-1/a))/(a+5)/(a+4)/(a+1)^6/(a+2)^3/(a+3)^2/(a+6))*s^7+O(s^8)

S:= [seq(subs(dsolve(sys union {op(v)}, numeric, output=listprocedure),[x(t),y(t)]),v=init)];

display([seq](pointplot(S(sqrt(2)*j/50),colour=[red,blue,gold,black]),j=1..50),insequence=true);

Since there are two independent variables (Y[1] and Y[2]), it's not a curve.  It's a mapping from a rectangle to the plane.

> U[1]:= ( -2 * ln(y1) )^(1/2) * cos(2*Pi* y2);
 U[2]:= ( -2 * ln(y1) )^(1/2) * sin(2*Pi* y2);
 plot([seq([U[1],U[2],y2=0 .. 1], y1=0.1*[$1..10]),
       seq([U[1],U[2],y1=0.1 .. 1], y2=0.1*[$1..10])],colour=[red$10,blue$10], scaling=constrained, axes=box);

If correlation of X(i) and X(i+1) is 0.9 and correlation of X(i+1) and X(i+2) is 0.9, then correlation of X(i) and X(i+2) is at least 2*(0.9)^2 - 1 = 0.62.  So you are requesting the mathematically impossible.

However, for a first-order autocorrelation of alpha and a second-order autocorrelation of beta, where  1 > alpha > -1, 1 > beta > 2 alpha^2 - 1, you could do this:

x[1]:= r[1]:
x[2]:= alpha*r[1] + sqrt(1-alpha^2)*r[2]:
for i from 3 to n do
  x[i] := (alpha^2-beta)/(alpha^2-1)*x[i-2] + alpha*(beta-1)/(alpha^2-1)*x[i-1] + 
      sqrt((beta-2*alpha^2+1)*(1-beta)/(1-alpha^2))*r[i]
end do:

Hint: multiply the matrix by a certain vector.

Maple does have trouble with integrals such as this.  See e.g. www.mapleprimes.com/forum/integratingcossintheta.

Since nobody else has tackled this, here's my attempt using the maplet elements MathMLEditor and MathMLViewer.

> with(Maplets[Elements]):

 maplet := Maplet([
  [Label("Enter a mathematical expression in the variable x")],
  [MathMLEditor[MMLE1](outputformat=content)],
  [Label("Enter a value for x")],
  [TextField[TF1](20)],
  [Button("Evaluate",Evaluate('MMLV1' = 'ff()'))],       
  [MathMLViewer[MMLV1]()],
  [Button("Done", Shutdown())]
    ]): 

ff := proc()
   local S, expr,v;
   S:= Maplets[Tools][Get](MMLE1);
   expr:= MathML[ImportContent](S);
   v:= parse(Maplets[Tools][Get](TF1));
   evalf(eval(expr,x=v));
  end proc:

Maplets[Display](maplet);

Note that by right-clicking in the top window you get access to palettes.

The standard way to interpret a^b^c in mathematics is a^(b^c), because (a^b)^c could simply be written as a^(b*c) (when this simplification is valid) .  Now 2^(56^89) is an enormous number (more than 10^155 decimal digits).  It's not at all surprising that Maple can't handle it.  By contrast, (2^56)^89 = 2^(56*89) = 2^4984 has only 1501 digits.