restart;
interface(version); 

fname:="G:/data/tmp/ftest.mpl":
f:=proc(n)
  local i,k;
  for i to n do f(i):=parse(cat("", seq(k mod 10, k=1..10*i))) od;
  f(n);
end:

f(80):
length(%);  # 800, ok
save f, fname;
f:='f':
#restart;
read fname:
f(80):
length(%);  # 100  ???

one way might be to use solve or PDEtools:-Solve

restart;
eq:=5*x^3-x^2+x-1=0;
PDEtools:-Solve(eq,x)

This is how I do these things.  

Rule of thumb, always use PDF files for images if possible. If not, export to .eps/ps, then convert the eps to PDF, then use that in the incluegraphics command in Latex. Latex compilers (pdflatex and lualatex, etc...) can read pdf image file for many years now.

I do not use the Maple worksheet GUI to export anything, as I do all this in code. Something like this

p:=plot(.....); #your plot command
plotsetup(ps, plotoutput="file_name_here",plotoptions = `noborder`);
print(p);
plotsetup(default):

The above will create file_name_here.ps. Choose any name you want for the file.

Now I run the commands

epspdf filename.ps
pdfcrop --margins 10  filename.pdf  filename.pdf

The first command convert the .ps file to pdf. The second crops it.

I've been doing this for years in all my Latex. It works very well for me. 

I just applied these to your image

with(plots):
scheme1 := ["zgradient",["Blue","Cyan","Green","Yellow","Orange","Red"]]:

P1:=plot3d(x*y, colorscheme=scheme1,style=surfacecontour ): 
P2:=plot3d(x*y,colorscheme=scheme1,style=point,symbol=asterisk): 
final_plot:=display({P1,P2}); #I want to export this figure to .eps
plotsetup(ps, plotoutput="file_name_here",plotoptions = `noborder`);
print(final_plot);
plotsetup(default):

Here is the result

\documentclass[11pt]{article}
\usepackage{graphicx}
\begin{document}
This is my graphics generated in Maple 2020.2 

\includegraphics[width=0.9\textwidth]{file_name_here}
\end{document}

In Mathematica, this is much easier, since Mathematica Export command can export plots and other Graphics to PDF directly.

May be in Maple 2021 exporting to PDF using the export() command will be supported.   I have no idea why Maple still does not support PDF format in the export command. PDF is the most universal format there is.

https://www.maplesoft.com/support/help/Maple/view.aspx?path=Export

foo5.pdf

There are all build in. But you could always make a function and make your own. For example for Hessian, something like

my_hessian:=proc(f::algebraic,vars::list(name))
 local n:=numelems(vars);
 local A:=Matrix(n,n);
 local i::posint,j::posint;
 for i from 1 to n do
   for j from 1 to n do 
       A[i,j]:=diff(diff(f,vars[i]),vars[j]);
   od;
od;
return A;
end proc;

You call it as

my_hessian(cos(x*y),[x,y])

Compare to

VectorCalculus:-Hessian(cos(x*y),[x,y])

It is ofcourse better to use the build in ones, as those have been tested over many years and have more options.

edit: changed function to algebraic. and vars also. Thanks to Carl Love reply below. 

It is better to make an attempt, and then ask for help when stuck. This way you'll learn better.

But I will make an attempt at part (a) for now.

f:=x->sin(x);
g:=x->cos(x);
plot([f(x),g(x)],x=0..Pi/2);

To find area bounded

h:=x->`if`(f(x)>g(x),f(x),g(x));
plot(h(x),x=0..Pi/2,view=[default, 0..1],filled=true);

To find the above area, use int. 

evalf(int(g(x),x=0..Pi/4)+int(f(x),x=Pi/4..Pi/2))

                1.414213562

One possible way

subs(exp(a)=p,expand(exp(-a)))

 Would you tell me what steps Maple might have gone through in order to come to the answer?

Maple does not show steps of solution. You might try the new Maple Learn application. (it can show step by step solution for some commands)

For Maple, you could look at trace by adding the command infolevel. For your problem it says

infolevel[solve]:=5;
eq_VR := VR = solve(eq_VS_m1, VR) assuming VS > 0, VR> 0;


solve: Warning: solve may be ignoring assumptions on the input variables.
solve:
             {VR::(RealRange(Open(0), infinity)), 


Main: Entering solver with 1 equation in 1 variable
Dispatch: dispatching to Rename handler
Dispatch: renaming conjugate(Sr) = _S000003
Recurse: recursively solving 1 equations and 1 inequations in 1 variables
Dispatch: handling polynomials of the form a*x^n-b
Dispatch: handling a single polynomial
Main: solving successful - now forming solutions
Main: Exiting solver returning 1 solution

Quadratic: solving a quadratic polynomial explicitly

So it looks like from the trace it used the standard quadratic formula.

I think the error  Error, (in PDEtools:-Library:-NormalizeBoundaryConditions)  is a side error, due to using real numbers with exact solver. Maple can't solve this PDE.

Rule of thumb. Do not use real numbers with exact solver. Rationalized all values.  After doing this, the error goes away, and pdsolve returns just () with no solution. Another thing, I do not like using eval to setup BC. Better using D[] directly.

Also, for the Laplacian in cylinderical coordinates, you can use VectorCalculus:-Laplacian to build it, instead of typing it yourself. Less change of error.

Also you had h:=(x,t)->0. But there is no x in this problem. I just set h:=0.

restart;
a := 35.5*0.0254: a:=convert(a,rational);
b := 36.5*0.0254: b:=convert(b,rational);
L := 0.0254*30*12: L:=convert(L,rational);
alpha := 0.000026972: alpha :=convert(alpha,rational);
Psi_s := 0.0440: Psi_s := convert(Psi_s,rational);
flux_b := 5.1/12: flux_b := convert(flux_b,rational);

h := 0:
the_laplacian_in_cylinderical :=VectorCalculus:-Laplacian( u(r,z,t), 'cylindrical'[r,phi,z] );
pde := diff(u(r, z, t), t) = alpha*the_laplacian_in_cylinderical + h;

bc:= u(r,L,t)=Psi_s, 
     -alpha*(D[2](u)(r,L,t))=0,
     -alpha*(D[1](u)(a,z,t))=0,
     -alpha*(D[1](u)(b,z,t))=flux_b;

ic:=  u(r, z, 0) = z*Psi_s/L;

#infolevel[pdsolve]:=5;
sol := pdsolve([pde,bc,ic],u(r, z, t)) ;

You might want to make sure your IC and BC are consistent. You IC says u=z*PSI/L. This means u_z = PSI/L, a constant. This is at t=0, but for all z and all r. But one of your BC says -apha*u_z = 0 (at z=L). Since alpha is not zero, this means u_z = 0 (at z=L). Which is conflicting with the initial condition. I do not know now if this is why the PDE was not solved.

Either way, Maple simply can't solve this analytically. No error generated.

Maple 2020.2 on windows

Your IC and BC are inconsistent.

You have 

ic := u(x, 0) = cos(x), u(0, t) = 0;

First condition says that for all x, at t=0, u=cos(x). 

The second says at x=0, for all t, u=0.

But at x=0, from the first conditions, it says u=1, since cos(0)=1 which is not zero according to second condition.

I could not understand your worksheet. But to solve this pde you could do the following

restart;
pde := diff(u(x,t),t) = k*diff(u(x,t),x$2)+h(x,t);
bc  := u(0,t) = 10,u(1,t) = 20;
ic  := u(x,0) = 60*x -50*x^2+10;
sol:=pdsolve([pde,bc,ic],u(x,t))

For specific h(x,t)

sol:=simplify(subs(h(x,tau)=0,sol))

you could see example in

https://rosettacode.org/wiki/Associative_array/Creation#Maple

"Maple tables are hashed arrays. A table can be constructed by using the table constructor."

T := table( [ (2,3) = 4, "foo" = 1, sin(x) = cos(x) ] );
          T := table(["foo" = 1, sin(x) = cos(x), (2, 3) = 4])
 
> T[2,3];
                                   4
 
> T[sin(x)];
                                 cos(x)
 
> T["foo"];
                                   1

if it for display, you could always do

a*z^``(-1)

I do not do plotting much in Maple, so this might not be the optimal solution. it seems to lose the minor ticks though.  There might be a way to also keep these there. I think the problem is that Maple needs a little bit more space, and so it added ticks by skipping every other value. that is why the odd values 5 and -5 do not show up.

R:=10:
r:=5:
plot3d( [ ( R+r(u,v)*cos(v))*sin(u),
          ( R+r(u,v)*cos(v))*cos(u),
            r*sin(v)
        ],
            u=0..2*Pi,
            v=0..2*Pi,style=patchnogrid,
            scaling=constrained,
            scaling=constrained,
coordinateview=[-15..15, -15..15,-5..5],tickmarks=[default,default, [seq(i,i=-5..5,2)]],
lightmodel=light3, viewpoint = circleleft);

Or you could increase the view to coordinateview=[-15..15, -15..15,-6..6] to give it more space? 

Maple uses last name evaluation for number of expressions. But this does not apply to things that are rtable based. This includes Matrix, Vector, and Array etc...

last name evaluations is what you are asking for. This means, when you write A:=C; and then assign some value to C, then when you use A again, you expect it to have same value you assigned to C. This is because Maple used last name evaluation to A.

But this does not work with rtables.  That is why when you print the matrix again, you did not see its value changed after you changed m

I do not know if it is possible to change this. But you can always use eval or subs

restart; 
u:=<-2+m,3+m>;
m:=5;
u;

type(u,'last_name_eval')

false

You can now use eval to force value of m to be used inside rtable

unassign('m');
u:=eval(u,m=5)

Look up last name evaluation in help for more information.

I get slighlty different result from that shown.

f:=t->862.979774+1825.011621/((1+1.54943476*e^(-677.0138344*t))^(131537/203808));
convert(f(t),rational)

Copied the integrand to Mathematica and it gives

Copied the result from Mathematica back to Maple

expr:=`(1450669*t)/
  1681 + (4459704101248*
    Hypergeometric2F1[131537/203808, 1, 335345/203808, 
          (36887*E^((342569*t)/506))/(57154 + 
        36887*E^((342569*t)/506))])/
     (1067737900495*(1 + 57154/(E^((342569*t)/506)*36887))^(131537/
       203808))`;
maple_expr:=MmaTranslator:-FromMma(expr);
lprint(maple_expr);


1450669/1681*t+4459704101248/1067737900495*hypergeom([131537/203808, 1],[335345
/203808],36887*exp(342569/506*t)/(57154+36887*exp(342569/506*t)))/(1+57154/
36887/exp(342569/506*t))^(131537/203808)

Compare to answer by Axel 

Here is plot on top of each others, which shows some difference. I am not sure which is the correct one

Axel:=2.695678478*hypergeom([1., 1., 1.645396648],[2., 2.],-1.549434760*exp(-677.0138340*t))*exp(-677.0138340*t)+.95671881+1887.991395*t;
plot([maple_expr,Axel],t=-0.01 .. 0.01)