because it relies so much on different files/worksheets on OP's machine.

I suppose it *might* be a version or interface issue.

The attached, using exclusively 1-D input, running Maple2015 on Win7 64-bit runs just fine for me.

NB in order to execute this worksheet properly, the path/filename entries will have to be changed to something appropriate for your machine in both the savelib() and libname() commands

unitsLibs.mw

If I understand correctly, then the attached should do what you want

plotNcomb.mw

To change the "number of atoms", just change the value at `Natoms:=`

Use rsolve(), on your "general" recursion relations

See the attached

recurSol.mw

restart;
a:=[1,2,3,4,5];
#
#Numerically
#
add(j*10^(numelems(a)-j), j in a);
#
#As string
#
cat(seq(convert(j, string), j in a));

Observation1: with a1=2, a2=2, a3=2, all three of the rational fractions are identically zero,  so the function value is 1

Observation2: for any other values of the parameters, each of the fractional terms is positive, so function is posVal1+posVal2+posVal3+1

Hence the minimum value of the function is 1, when a1=2, a2=2, a3=2

Observation3: for "large values of all three parameters, each of the fractional terms has a numerator which is quartic in 'a', and each denominator is quadratic in 'a', so each fractional term is quadratic in 'a'.

Hence the maximum value of the function is infinity when a1=infinity, a2=infinity, a3=infinity

The condition i+j != 0 on the inner summation seems to be causing some kind of issue.

i=0, j=1 is a valid combination of indices: contributing a value of 1 to the final sum
i=1, j=0 is also a valid combination of indices: contributing another value of 1 to the final sum.

Hence the sum has to be greater than 2: any answer less than 2, is therefore wrong.

I failed to get any kind of "closed-form" solution to this problem, and had to resort to the use of add(), rather than sum(). After a bit of experimentation, setting 1000 as the upper value on the indices, seemed a pretty good balance between accuracy and execution time. My code for this problem is

  restart:
  t1:=time():
  ulim:= 1000:
  p:= Array(0..ulim, i->evalf(i^4)):
  add
  ( add
    ( 1/(p[i]+p[j]),
      j =`if`
           ( i=0,1,0 )..ulim
    ),
    i=0..ulim
  );
  t2:=time()-t1;

which return the answer 2.926377298 in about 1.25secs on my machine.

Precalculating all the relevant fourth powers (as in the above), and using lookups is about 2X faster (for ulim=1000) than using the really simple approach below

  restart:
  t1:=time():
  ulim:=1000:
  add(add(1.0/(i^4+j^4), j=`if`(i=0,1,0)..ulim), i=0..ulim);
  t2:=time()-t1;

1. Funny how questions involving units don't seem to get many answers?
2. I personally avoid the use of any "Units" package in almost any software (including Maple). I just make sure that all of my quantities are expressed in MKS (I'm European, so it is natural), then I know what the "units" of any answer will be
3. Having got that of my chest I cannot actually see anything wrong with your worksheet - it really *ought* to function, without an irrelevant (and incorrect) complaint about units
4. To illustrate the statement in (3) above, I commented out the last couple of lines in your original worksheet. I replaced them with a "step-by-step" way to get the answer for the quantity you want. This appears to function, provided that the expected units of the quantity f_co is a "pressure", so that units of N/m² are correct. See the attached file for how I achieved this
5. If the units of the final answer are correct (and the numerical value is OK) then I think the problem is with Maple, not with your original code

UnitProb.mw

You can use the plottools[exportplot]() command

Interpret the "integral" as a simple Riemann Sum (see the Wikipedia article)

rho:=[0.045, 0.0459, 0.0564, 0.05689, 0.06015, 0.06235, 0.0654, 0.0687, 0.07012, 0.07251];
eta__k:=[1.15, 1.256, 1.56, 1.85, 1.86, 2.01, 2.35, 2.56, 2.86, 2.901];
f:=x->rho[x]*x^4*eta__k[x]^2;
leftRiemannSum:=add(f(x),x=1..9);
rightRiemannSum:=add(f(x),x=2..10);
TrapezoidalRule:=(1/2)*( 2*add(f(x),x=1..10)-f(1)-f(10) );

will return

leftRiemannSum := 6993.111707
rightRiemannSum := 13095.34950
TrapezoidalRule := 10044.23060

Since your function is montonically increasing, one would expect the left Riemann sum to be an underestimate, the right Riemann sum to be an overestimate and the trapezoidal rule to be "pretty close".

Note that the result from the trapezoidalRule is close (within 2%) to that obtained by Kitonum

Have you considered

Student[Calculus1][CriticalPoints]()???

whic would seem to cover all of your requirements

I'm too old to work with fancy features like 2D Math Input, equation references etc, but you could try the following. In 1-D inpu and no equation references, it just works

restart;
with(PDETools):
declare( F(r,s), H(r), G(s) ):
sys2 := [ -diff(F(r,s),r,r) + diff(F(r,s),s,s) + diff(H(r),r) + diff(G(s),s) + s = 0,
              diff(F(r,s),r,r) + 2*diff(F(r,s),r,s) + diff(F(r,s),s,s) - diff(H(r),r) + diff(G(s),s)-r = 0
           ];
sol:=pdsolve(sys2);

So if I wanted 'x' in the first column, sin(x) in the second column and x^2 in the third column, I would probably use something like the following

    restart;
#
# Generate the data as a listlist
#
    A:=[seq( [j, sin(j), j^2], j=-5..5, 0.1)]:
#
# Make sure everything has evaluated to a float
#
    A:=map(evalf, A);
#
# define an appropriate filename. OP will have to ensure that this
# path, etc is appropriate for  their machine - so changes to the
# pathname will be necessary
#
   fname:="J:/Users/TomLeslie/myMaple/writeTest.dat":
#
# Actually write the data
#
   writedata(fname, A);

(1) How do i make output (that blue part i get after pressing Enter) go in one line?

Pretty sure the above can't be done.

(2) And also - is there a way to hide output?

If you don't wnat the output of a statement to display, terminate the statement with a colon, rather than a semicolon.

A:=[2,4,6,8,10,12];

ListTools[Search](6,A);

will return the first index with 6 as value

If you want all occurences, for example in the list

A:=[2,4,6,8,10,6,12,6], then

ListTools[SearchAll](6,A);

will return a sequence of all indices which have value 6

The expression for your cusp catastrophe, ie x^4+x^2*z+x*y, produces zero whenever x=0. This is why you have the "horizontal plane" in your output plot - Maple is plotting what you asked for!

I have seen several different equations which define the cusp catastrophe and avoid this problem, but if you want to stick with your definition, then just avoid x=0: for example you could try

p1:=implicitplot3d( x^4+x^2*z+x*y,
                            x = -7 .. -1e-06,
                            y = -50 .. 50,
                            z = -20 .. 20,
                            axes = boxed,
                            style = surfacecontour,
                            grid = [60, 60, 60],
                            orientation = [180, 8, 175],
                            'transparency' = .1
                        ):
p2:=implicitplot3d( x^4+x^2*z+x*y,
                            x = 1e-06 .. 7,
                            y = -50 .. 50,
                            z = -20 .. 20,
                            axes = boxed,
                            style = surfacecontour,
                            grid = [60, 60, 60],
                            orientation = [180, 8, 175],
                            'transparency' = .1
                         ):
plots[display]([p1,p2]);

There are many other ways to do this, but probably not with your starting equation, and its x=0 issue