> _EnvAllSolutions:= true;
   solve(sin(x)=1/2);

1/6*Pi+2/3*Pi*_B1+2*Pi*_Z1

Note that _B1 is a "binary" variable, i.e. 0 or 1, and _Z1 is an integer variable.

> solve(cos(4*x)=sin(x));

1/2*Pi+2*Pi*_Z2, 1/2*Pi+2*Pi*_Z3, -1/6*Pi+2*Pi*_Z4,

-5/6*Pi+2*Pi*_Z4, arctan(4*(1/4*5^(1/2)-1/4)/(10+2*5^(1/2))^(1/2))+2*Pi*_Z5,

-arctan(4*(1/4*5^(1/2)-1/4)/(10+2*5^(1/2))^(1/2))+Pi+2*Pi*_Z5, arctan(4*(-1/4-1/4*5^(1/2))/(10-2*5^(1/2))^(1/2))+2*Pi*_Z6, -arctan(4*(-1/4-1/4*5^(1/2))/(10-2*5^(1/2))^(1/2))-Pi+2*Pi*_Z6

What precisely do you mean by a Gaussian fit to a histogram?  Do you mean that your data are a random sample from some unknown distribution, and you want a normal distribution that approximates this?  The usual practice would be to use the sample mean and sample standard deviation (the one with N-1 instead of N).  This is provided by StandardDeviation in the Statistics package.  But it's a Gaussian fit to the data, not a Gaussian fit to a histogram. Do you then want to show the corresponding Gaussian curve and the histogram in the same plot?

By the way, you should be using Statistics rather than the old stats package.

Thus:

> with(Statistics):
   dataset := Vector([7.6, 8.1, 4, 1, 2.1, 2.1, 2.1, 3.3, 5, 4, 3, 3, 2.5, 2.5, 2.9,
      2.9, 2.9, 1.5, 1.75], datatype = float);
   Q := Histogram(dataset, maxbins = 10, binwidth = .5, frequencyscale = absolute, 
      range = 0 .. 10, color = grey);
   mu:= Mean(dataset);
   sigma:= StandardDeviation(dataset);
   G:= plot(Count(dataset)*0.5*PDF(Normal(mu,sigma),t), t=0..10):
   plots[display]([Q,G],  title = "Histogram and Gaussian", axes = boxed, 
      labels = ["Droplet diameter (nm)", "Number of droplets"], 
      labeldirections = [horizontal, vertical]);

As you can see, it doesn't look like a very good fit, mainly I think because of those outliers at 7.6 and 8.1.

Well, your first error seems to be that the procedure A_y is calling  ArrayInterpolation(Z, Ay, Ratio) with Ratio not a numerical value but a number with units attached.  ArrayInterpolation doesn't know anything about units: it needs numerical values.

> series(a/x + b/x^2, x);

Use semilogplot from the plots package.  For example:

> plots[semilogplot](sin(ln(x)),x=1e-7 .. 10);

Automatic simplification distributes rational numbers over sums.  If you want something that looks like (a[2]+a[3])/2, you can try (in Standard GUI)

> (a[1] + a[2])/`#mn("2")`;

However, that's really a name rather than an integer in the denominator, so you may have trouble if you try to use this in further computations.

You might try

> kernelopts(opaquemodules=false);
  showstat(ModularSolver:-Methods:-Mod2:-Solve);

I don't believe it is possible to do any of these, or anything that is not listed in the help page as either accessible by DocumentTools[SetProperty] or on the Component Properties menu.

> with(plots):
  implicitplot([x-exp(-y^2)=0, sin(5*x) + 2 - (y-1)^2 = 0], x = -2 .. 2, y = -2 .. 2, colour=[red,blue], gridrefine=3);

Try using the range option.  Thus compare:
 

> Data:= [[1,2,3,4,5,6,7,8,9]];
   listdensityplot(Data, colorstyle=HUE);
   listdensityplot(Data, colorstyle=HUE,range=0..15);  
   

A:= sort(Array(1..3, i->M[i]), (a,b) -> (a[1]<=b[1]));
Matrix(3,2,(i,j)->A[i][j]);

The command is odeplot, not odeplots. 

If the system is third order or higher in X[5], dsol should provide values of that second derivative, and then odeplot
will be able to plot it.  If not, you'd need to find an expression for the second derivative in terms of variables and their first derivatives. You can use rifsimp in the DEtools package to help with this.  For example, consider this system:

> sys:= {D(x)(t)=y(t)^2 + x(t), (D@@2)(y)(t)=x(t)};
    ics:= {x(0) = 1, y(0) = 0, D(y)(0) = 0};
    
> DEtools[rifsimp](convert(sys union {(D@@2)(x)(t)=w(t)},diff),[[w],[x,y]]); 

TABLE([Solved = [w(t) = 2*y(t)*diff(y(t),t)+y(t)^2+x(t), diff(y(t),`$`(t,2)) = x(t), diff(x(t),t) = y(t)^2+x(t)]])

> d2x:= subs(%[Solved], w(t));

d2x := 2*y(t)*diff(y(t),t)+y(t)^2+x(t)

> plots[odeplot](dsol, [t, d2x], t = 0 .. 1);

There is no method called "implicit".  If you're talking about numerical methods, the rosenbrock method is described as an "Implicit Rosenbrock third-fourth order Runge-Kutta method" and the mebdfi method is a "Modified Extended Backward Differentiation Equation Implicit method".  Also, most of the numerical methods support the option implicit = true, which means (for a system linear in the leading derivatives) don't solve symbolically for those leading derivatives, but instead solve numerically when they are needed.

Or are you referring to the implicit option when solving symbolically, which means avoid trying to get the solution in explicit form?

One way (which would be OK in small examples such as this) would be to generate a sequence of random matrices and accept the first one where the condition is satisfied.  Thus:

>  with(LinearAlgebra):
    do
      A:= RandomMatrix(3,2,generator=1..3);
      if nops(convert(Column(A,1), set)) = 3 then break end if;
      end do:
   A;
                               [3    1]
                               [      ]
                               [2    1]
                               [      ]
                               [1    1]

Another method would be to use randcomb and randperm in the combinat package to generate the first column of A.  This might be better in cases where very few of the random matrices would satisfy the condition.  Thus for a 50 x 2 matrix with entries from 1 to 100 and no repeated elements in the first column:
 

> with(combinat):
    A:= < convert(randperm(randcomb(100,50)),Vector) |
            RandomVector(50, generator=1..100) >;

1) You have an extra right parenthesis at the end of the x[bar] computation, so you get "unable to match delimiters".

2) y is a Vector; there are elements y[1] to y[8], but there is no such thing as y[bar].  This would also affect x[bar] if you fixed the error from (1).  You could convert these to atomic identifiers (as in (3)), or use different names such as xbar and ybar.

3) You gave values to what look like M[x] and M[y], but are really symbols `#msub(mi("M"),mi("x"))` and  `#msub(mi("M"),mi("y")).
But then in the Q[i] computation you use M[x] and M[y].  You need to highlight each of those, right-click and choose

2D Math, Convert to, Atomic Identifier