Start with the following command:

> plots[setoptions](axesfont = [Courier, 20]);

... and then call up the Curve Fitting Assistant; it should use 20-point Courier for the numbers on the axes.

You see only one plot because p1, p2 and p3 are the same.  Only the names of the variables are changed.
Maple doesn't automatically put x, y and z on the first, second and third axes, it takes them in the order that you
give them.

To get my plot, all you should need is this:

 

> cylx := y^2+z^2-1:
  cyly := z^2+x^2-1:
  cylz := x^2+y^2-1:
  expr := max(cylx, cyly, cylz);
  plots:-implicitplot3d(expr, x=-1..1, y=-1..1, z=-1..1, numpoints=10\000);

You see only one plot because p1, p2 and p3 are the same.  Only the names of the variables are changed.
Maple doesn't automatically put x, y and z on the first, second and third axes, it takes them in the order that you
give them.

To get my plot, all you should need is this:

 

> cylx := y^2+z^2-1:
  cyly := z^2+x^2-1:
  cylz := x^2+y^2-1:
  expr := max(cylx, cyly, cylz);
  plots:-implicitplot3d(expr, x=-1..1, y=-1..1, z=-1..1, numpoints=10\000);

> Results := fsolve(...);
  result_liq := eval(psi_liq, Results);

In that case you can use writeto or appendto.  For example:

> interface(echo=0);
  writeto("c:/temp/myfile.txt");
  f:= (x^2+a^2)/(x^2+b^2) ;
  writeto(terminal);

and the file will now contain

                                       2    2
                                      x  + a
                                 f := -------
                                       2    2
                                      x  + b

In that case you can use writeto or appendto.  For example:

> interface(echo=0);
  writeto("c:/temp/myfile.txt");
  f:= (x^2+a^2)/(x^2+b^2) ;
  writeto(terminal);

and the file will now contain

                                       2    2
                                      x  + a
                                 f := -------
                                       2    2
                                      x  + b

> f:= x -> exp(3*x)/((1+exp(x))^5);
  ff:= int(f(x)*f(y-x),x=-infinity..infinity) assuming y::real;


  ff := 1/6 exp(3 y) (-160 exp(3 y) - 25 exp(4 y) + 160 exp(y)

         + 96 y exp(3 y) + 6 y exp(4 y) + 216 y exp(2 y)

         + 96 y exp(y) + 25 + 6 y)/(exp(9 y) - 9 exp(8 y)

         + 36 exp(7 y) - 84 exp(6 y) + 126 exp(5 y) - 126 exp(4 y)

         + 84 exp(3 y) - 36 exp(2 y) + 9 exp(y) - 1)

> series(ff,y=0,20);

                             2            3     71     4     137    5
  1/630 + 1/1260 y - 1/3960 y  - 2/10395 y  + ------- y  + ------- y
                                              4324320      5405400

                        6      547     7      2039      8
         - 1/129729600 y  - --------- y  - ----------- y  +
                            227026800      17643225600

           4679      9       290327      10      11
        ----------- y  + -------------- y   + O(y  )
        25261891200      19199037312000

> Digits:= 100; # a bit of overkill
  ff1:=numapprox[minimax](ff,y=-0.2 .. 0.2, 8, 1, 'maxe');
  ff2:= evalf[15](ff1);

  ff2 := 0.00158730158730159 + (0.000793650793634240 + (

        -0.000252525252531389 + (-0.000192400186878993 + (

        0.0000164187674671730 + (0.0000253445280276351 + (

                             -8
        -0.777281641667611 10

                                 -5                       -6
         + (-0.239280581666087 10   - 0.113915499760773 10   y) y) y)

        y) y) y) y) y

  
 

>  evalf[15](maxe);

.367783513840569e-15

So ff2 is a good approximation to ff for -0.2 < y < 0.2, with maximum error less than 0.4 * 10^(-15).
Outside that interval, hopefully, the exact solution is not too bad numerically.

> f:= x -> exp(3*x)/((1+exp(x))^5);
  ff:= int(f(x)*f(y-x),x=-infinity..infinity) assuming y::real;


  ff := 1/6 exp(3 y) (-160 exp(3 y) - 25 exp(4 y) + 160 exp(y)

         + 96 y exp(3 y) + 6 y exp(4 y) + 216 y exp(2 y)

         + 96 y exp(y) + 25 + 6 y)/(exp(9 y) - 9 exp(8 y)

         + 36 exp(7 y) - 84 exp(6 y) + 126 exp(5 y) - 126 exp(4 y)

         + 84 exp(3 y) - 36 exp(2 y) + 9 exp(y) - 1)

> series(ff,y=0,20);

                             2            3     71     4     137    5
  1/630 + 1/1260 y - 1/3960 y  - 2/10395 y  + ------- y  + ------- y
                                              4324320      5405400

                        6      547     7      2039      8
         - 1/129729600 y  - --------- y  - ----------- y  +
                            227026800      17643225600

           4679      9       290327      10      11
        ----------- y  + -------------- y   + O(y  )
        25261891200      19199037312000

> Digits:= 100; # a bit of overkill
  ff1:=numapprox[minimax](ff,y=-0.2 .. 0.2, 8, 1, 'maxe');
  ff2:= evalf[15](ff1);

  ff2 := 0.00158730158730159 + (0.000793650793634240 + (

        -0.000252525252531389 + (-0.000192400186878993 + (

        0.0000164187674671730 + (0.0000253445280276351 + (

                             -8
        -0.777281641667611 10

                                 -5                       -6
         + (-0.239280581666087 10   - 0.113915499760773 10   y) y) y)

        y) y) y) y) y

  
 

>  evalf[15](maxe);

.367783513840569e-15

So ff2 is a good approximation to ff for -0.2 < y < 0.2, with maximum error less than 0.4 * 10^(-15).
Outside that interval, hopefully, the exact solution is not too bad numerically.

For just the intersection:

expr := max(cylx, cyly, cylz);

plots:-implicitplot3d(expr, x=-1..1, y=-1..1, z=-1..1, numpoints=10\000);

For just the intersection:

expr := max(cylx, cyly, cylz);

plots:-implicitplot3d(expr, x=-1..1, y=-1..1, z=-1..1, numpoints=10\000);

Why not Matrices?  It seems to work.

> A:= Matrix(1,1):
  A(3,6):= 7:
  A;

                     [0    0    0    0    0    0]
                     [                          ]
                     [0    0    0    0    0    0]
                     [                          ]
                     [0    0    0    0    0    7]

Why not Matrices?  It seems to work.

> A:= Matrix(1,1):
  A(3,6):= 7:
  A;

                     [0    0    0    0    0    0]
                     [                          ]
                     [0    0    0    0    0    0]
                     [                          ]
                     [0    0    0    0    0    7]

Sorry, I should have tried this out first.  The problem is that the for loop leaves i with a value (which wouldn't be a problem if I used seq, but for some reason I felt like using $).

plot([h[j] $ j=1..10], x = 0 .. 2*Pi);

would work.

Sorry, I should have tried this out first.  The problem is that the for loop leaves i with a value (which wouldn't be a problem if I used seq, but for some reason I felt like using $).

plot([h[j] $ j=1..10], x = 0 .. 2*Pi);

would work.

with(Statistics):
Phi:= unapply(CDF(Normal(0,1),x), x);

Phi := x -> 1/2+1/2*erf(1/2*x*2^(1/2))