restart:

  ra:=2: b1:=1.41: na:=0.7: we:=0.5: eta[1]:=4*0.1: d:=0.5/1:
  xi:=0.1: m:=na: ea:=0.5: pr:=21: gr:=0.1: R:=0.9323556933:

  PDE1:=ra*(diff(f(x,t),t))=+b1*(1+ea*cos(t))+(1/(R^2))*((diff(f(x,t),x,x))+(1/x)*diff(f(x,t),x));
  IBC:= {D[1](f)(0,t)=0,f(1,t)=0,f(x,0)=0};

#
# Solve the PDE then use the returned methods 'plot', 'plot3d',
# and 'animate' to produce various plots
#
  sol :=  pdsolve({PDE1}, IBC, numeric) :
  sol:-plot(f(x, t), t = 1.2, linestyle = "solid", title = "Velocity Profile", labels = ["r", "f"]);
  sol:-plot3d(f(x, t), x=0..1, t=0..2, style=surface, color=cyan );
  sol:-animate( f(x,t), x=0..1, t=0..2);

#
# Use the 'value' method to facilitate computation of
# assorted numerical values
#
# Check which quantities are "immediately available"
#
  sol:-value(f(x,t), output=listprocedure);
#
# Aow get some numerical values from this basic method
#
  sol:-value(f(x,t))(0.5,0.5);
#
# For ease of use, one can assign the 'f(x,t)' procedure to,
# the name 'fN' and then compute the values, as in
#
  fN:=eval( f(x,t), sol:-value(f(x,t), output=listprocedure)):
  fN(0.5,0.5);

#
# The solution module returned by pdsolve() does not
# contain any derivatives, so these have to be computed
# explicitly. The simplest method is to use the 'D'
# operator.
#
# Differentiation wrt 'x' evaluated at x=0.2, t=1.2
#
  D[1](fN)(0.2, 1.2);
#
# Differentiation twice wrt 'x' and evaluate at x=0.2,
# t=1.2
#
  D[1,1](fN)(0.2, 1.2);
#
# Differentiation wrt 't' evaluated at x=0.2, t=1.2
#
  D[2](fN)(0.2, 1.2);
#
# Plot the first and second derivatives of f(x,t) wrt 'x' for t=1.2
# Note the "glitch" in the second derivative
#
  plot( [ D[1](fN)(x, 1.2),
          D[1, 1](fN)(x, 1.2)
        ],
        x=0..1,
        color=[red, blue]
      );
#
# Plot the first and second derivatives of f(x,t) wrt 't' for x=0.5
#
  plot( [ D[2](fN)(0.5, t),
          D[2, 2](fN)(0.5, t)
        ],
        t=0..2,
        color=[red, blue],
        axes=boxed
      );
  

  M:= Matrix([ [ "x", "f(x,t)", "diff(f(x,t),x)", "diff(f(x,t),x,x)"],
                  seq( [j, fN(j, 1.2), D[1](fN)(j,1.2), D[1,1](fN)(j,1.2)], j=0.1..0.9, 0.1)
               ]
            );
 # ExcelTools:-Export( M, "C:/Users/TomLeslie/Desktop/pdeDat.xlsx")