The order of elements in a set is arbitrary.  Unfortunately, fsolve does not accept lists rather than sets (as solve does).

As for the second question:

S:= fsolve(...):

aHat:= subs(S,a);

bHat:= subs(S,b);

Searching for 1904128 in the On-line Encyclopedia of Integer Sequences yields

<http://www.research.att.com/~njas/sequences/A012696>

which is just your sequence, but no other information about it is listed except the definition (which is closely related to yours).

 But I might note that your numbers are also the derivatives of arctan(cos(2*x)) at x=0,

i.e. arctan(cos(2*x)) = arctan(tan(x+Pi/2)^2) - Pi/4

.

select(t -> (t mod 12 = 4), [$-100..100]);

What happens when xi is very large and negative?

>   convert(x*BesselI(0,x), FormalPowerSeries, x);

Sum(4^(-k)/k!^2*x^(2*k+1),k = 0 .. infinity)

>  Q:= eval(op(1,%),x=sin(phi));

Q := 4^(-k)/k!^2*sin(phi)^(2*k+1)

> J:= unapply(Int(sin(phi)^(2*k+1),phi=0..Pi),k);

Now get a reduction formula for J(k).

> with(IntegrationTools):
   Parts(J(k),sin(phi)^(2*k)) assuming k::posint;

-Int(-2*cos(phi)^2*sin(phi)^(2*k)*k/sin(phi),phi = 0 .. Pi)

> eval(%,cos(phi)^2=1-sin(phi)^2);

-Int(-2*(1-sin(phi)^2)*sin(phi)^(2*k)*k/sin(phi),phi = 0 .. Pi)

> V:= map(combine,Expand(%));

So the reduction formula is J(k) = 2*k*J(k-1) - 2*k*J(k).  For k = 0 we have

> value(J(0));

2

Solve the reduction formula with this initial condition:

> rsolve({R(k) = 2*k*R(k-1) - 2*k*R(k), R(0) = 2}, R(k));

1/GAMMA(k+3/2)*Pi^(1/2)*GAMMA(k+1)

So the term for k in the integral is

> simplify(eval(Q,sin(phi)^(2*k+1)=%)) assuming k::posint;

2/GAMMA(2*k+2)

Summing over all k:

> Int(sin(phi)*BesseI(0,sin(phi)),phi=0..Pi) = sum(%,k=0..infinity);

Int(sin(phi)*BesselI(0,sin(phi)),phi=0..Pi) = 2*sinh(1)

The "." operator used for Matrix multiplication is non-commutative and associative.

> a . (b . c) - (a . b) . c;

              0

However, I don't know how to implement the distributive laws.

You want a sample from a truncated normal distribution.  Now unfortunately truncated normal is not one of the already-defined distributions, but we can get around that.

Suppose you're interested in the normal distribution with mean mu and standard deviation sigma, truncated below at a and above by b (in your case, I think,
a = max(0, mu - 3*sigma) and b = mu + 3*sigma.  If FN is the CDF of the Normal(mu,sigma) distribution, then the CDF of the truncated distribution is

F(x) = (FN(x) - FN(a))/(FN(b)-FN(a))

Given a random variable U with uniform distribution in [0,1], we get one with this truncated normal distribution by solving F(X) = U, a <= X <= b.  Thus (choosing some numerical values), here's one way to get a sample of size 1000.
 

> with(Statistics):
  mu:= 5.0; sigma:= 2.0; 
  a:= max(0,mu-3*sigma); b:= mu+3*sigma;
  FN:= unapply(CDF(Normal(mu,sigma),x),x);
  F:= unapply((FN(x) - FN(a))/(FN(b)-FN(a)), x);
  U:= RandomVariable(Uniform(0,1));
  S:= map(u -> fsolve(F(x)=u, x=a..b), Sample(U, 1000));

You posted the same question three times, but didn't say what is S and what is T.   To see whether these sums are the same, consider (without using Maple) what are the terms for n= 1, 2, 3...

Your X and Y have different numbers of elements.  I'll take

> X := Vector([1, 2, 3, 4, 5], datatype=float):
Y := Vector([2, 3, 4.8, 10.2, 15.6, 30.9], datatype=float):

> F:= subs(dsolve({diff(y(x),x) + u*y(x)=v, y(0)=0}, y(x), 
     numeric, parameters=[u,v], output=listprocedure),y(x));
  P:= [seq(subs(XX=X[i], YY=Y[i], 
          proc(u,v) F(parameters=[u,v]); YY - F(XX) end proc),
        i=1..5)];
> Optimization:-LSSolve(P);


  (Z)..Warning :, Must check assignment of Global var. for ,

        proc(x)  ...  end proc


  (Z)..Warning :, Must check assignment of Global var. for ,

        proc(x)  ...  end proc


  (Z)..Warning :, Must check assignment of Global var. for ,

        proc(x)  ...  end proc


  (Z)..Warning :, Must check assignment of Global var. for ,

        proc(x)  ...  end proc


  (Z)..Warning :, Must check assignment of Global var. for ,

        proc(x)  ...  end proc


                                 [-0.408540690360691240]
             [0.836622896559720, [                     ]]
                                 [0.956095638634626100 ]

So the optimal parameter values seem to be u=-0.408540690360691240, v = 0.956095638634626100.

The warnings seem to be harmless; I would have expected interface(warnlevel=0) to remove them, but it doesn't.  I'm submitting that as an SCR.

For example, one of your expressions is

> diff(Wrbajas(t90),t90); 

(sum(A[i]*(ln(t90/K+a)/a)^i*i/(K*(t90/K+a)*ln(t90/K+a)), i = 1 .. n))*exp(A0+sum(A[i]*(ln(t90/K+a)/a)^i, i = 1 .. n))

> simplify(%) assuming positive;

(sum(A[i]*a^(-i)*(ln(t90+a*K)-ln(K))^(i-1)*i, i = 1 .. n))*exp(A0+sum(A[i]*a^(-i)*(ln(t90+a*K)-ln(K))^i, i = 1 .. n))/(t90+a*K)

> combine(%,ln) assuming positive;

(sum(A[i]*a^(-i)*ln((t90+a*K)/K)^(i-1)*i, i = 1 .. n))*exp(A0+sum(A[i]*a^(-i)*ln((t90+a*K)/K)^i, i = 1 .. n))/(t90+a*K)

 [ why is the Maple tag not working? ]

If you just plug Eq1 into Eq2, you get a differential equation for W(t) that still has R(t) in it as an unknown function.  The solution will be an integral involving that unknown function.  It's correct, but it may not be very useful if you don't know what R(t) is.  Instead, you should use dsolve on the system of two equations together, and Maple will tell you what both W(t) and R(t) are.  Of course, there is still an unknown function H(t), but I suppose you can deal with that one yourself. 

It's hard to help when we don't know what the question is.  What are you trying to do?

I'm confident it can be done, but I think you'll have to tell in more detail what you want the program to do.  Do you want something like this?

> with(Student[Calculus1]):
    ApproximateInt(x^2, x=0..1, partition=
     [seq((3/4)^(10-j),j=1..10)], method=left,output=plot);

I don't know what point you want to zoom in on for the first plot, but the way to do it is to specify new intervals for x,y, and z.

> ?plots,implicitplot