sol0:=dsolve( eval(eqn1,[M = 0.4, pr = 0.0441]),numeric):
Mvals:=[0.00, 0.02, 0.04, 0.06]:
prVals:=[0.044, 0.71, 2.0, 7.0]:
printf( "\n\n%s     %s     pr=%5.3f pr=%5.3f pr=%5.3f pr=%5.3f\n\n",
        "M", "f'(0)", prVals[1], prVals[2], prVals[3], prVals[4]
    );
for i from 1 by 1 to 4 do
    ans:=[];
    for j from 1 by 1 to 4 do
        sol:=dsolve( eval(eqn1,[M = Mvals[i], pr = prVals[j]]),
                     approxsoln=sol0, numeric, maxmesh=1024
                   ):
        ans:=[ans[],rhs~(sol(0)[4..5])];
    od:
    printf( "%4.2f  %7.5f   %7.5f  %7.5f  %7.5f  %7.5f\n",
            Mvals[i],ans[1][1], ans[1,2], ans[2][2], ans[3][2], ans[4][2]
          )
od:

M     f'(0)     pr=0.044 pr=0.710 pr=2.000 pr=7.000

0.00  0.33206   0.20126  0.84174  1.40700  2.52723

0.02  0.35884   0.46479  1.22516  1.84743  3.06165

0.04  0.38397   0.72850  1.60825  2.28665  3.59244

0.06  0.40771   0.99235  1.99105  2.72481  4.12010

  restart;
#
# Increase the calculation precision to
# see if that helps - it doesn't, really
#
  Digits:=20:
  with(LinearAlgebra):
  a:= Vector([1,  2,     3]):
  b:= Vector([1, -2,     3]):
  c:= Vector([3,  2,     1]):
  d:= Vector([5,  1,    -3]):
  e:= Vector([0,  1.01, -3]):
  GramSchmidt([a, b, c, d, e]);
  GramSchmidt([b, c, d, e, a]);
#
# Set the calculation precison back to
# the default, and then use "exact
# rationals" rather than floats
#
  Digits:=50:
  f:= x-> convert~(x, rational):
  GramSchmidt(f~([a, b, c, d, e]));
  GramSchmidt(f~([b, c, d, e, a]));

#
# Initialize, load the CurveFitting() and
# plots() packages
#
  restart;
  with(CurveFitting):
  with(plots):
#
# Define the points
#
  L:=[[2, 2], [12, 6], [37, 42], [49, 21], [73, 49], [91, 2]]:
#
# Fit a polynomial in x to these points
# (and if you don't like it exact rationals,
# then evaluate it in floats
#
  poly:=PolynomialInterpolation( L, x);
  evalf(poly);
#
# Plot the original datapoints
#
  p1:=pointplot( L,
                 symbol=solidcircle,
                 symbolsize=20,
                 title= "Points to Fit",
                 titlefont=[times, bold, 20]
               );
#
# Plot the fittted curve from the minimum
# to the maximum x-values define by the
# list of points
#
  p2:=plot( poly,
            x= min
               ( [ seq
                   ( L[j,1],
                     j =1..numelems(L)
                   )
                 ]
               )..
               max
               ( [ seq
                   ( L[j,1],
                     j =1..numelems(L)
                   )
                 ]
               ),
            title= "Fitted Curve",
            titlefont=[times, bold, 20]
          );
#
# Display the original points and the fitted curve
# on the same graph
#
  display( [p1,p2],
           title= "Points and Fitted Curve",
            titlefont=[times, bold, 20]
         );

  restart;
  with(ScientificConstants):
  GetElement( C );
  evalf( Element( C, density, units ) );

#
# Change the density of 'carbon' so that
# subsequent call will return the value
# appropriate for diamond, rather than
# graphite
#
  ModifyElement(C, density=[ value = 3.5,
                             uncertainty = undefined,
                             units = g/cm^3
                           ]
               );
  evalf( Element( C, density, units ) );

#
# Alternatively define the value 'field' of
# density to be a procedure which accepts either
# 'graphite' or 'diamond' and returns the
# appropriate value
#
  ModifyElement( C, names={graphite, diamond},
                    density=[ value = proc()
                                          if   args[1]='graphite'
                                          then return 2.2
                                          elif args[1]='diamond'
                                          then return 3.5
                                          fi;
                                      end proc ,
                              uncertainty = undefined,
                              units = g/cm^3
                          ]
             );
  evalf( Element( C, density(graphite), units ) );
  evalf( Element( C, density(diamond), units ) );

restart;
#
# Converted the following to assignments
# and replaced 'pi' with 'Pi'
#
  a:= 0.11e-1:           b:= 0.25e-2: t:= 0.25e-2:
  theeta:= (1/180)*Pi:   L:= 0.17e-1: R:= L/theeta:
  shearmodulus:= 314000:
#
# eq0 doesn't seem to be used?? Comment out
#
# eq0 := Ma = int
#             ( int
#               ( (x*sin(O)+b)*x^2*(cos(O)*cos(O)),
#                 O = 0 .. (1/2)*Pi
#               ),
#               x = 0.9e-2 .. 0.11e-1
#             )-9;
#
# Don't need the assign() statement, since I
# converted the parameter definitions above
# from equations to assignments. So comment
# out
#
#  assign( ('t', 'b', 'L', 'shearmodulus', 'theeta', 'a')
#           =
#          (0.25e-2, 0.25e-2, 0.17e-1, 314000, (1/180)*Pi, 0.11e-1)
#        );
#
#
# eq4 doesn't seem to be used anywhere, so comment
# out
#
#  eq4 := Mo1 = int
#               ( int
#                 ( 2*shearmodulus*(1+Beeta*theeta/L-1/(1+Beeta*theeta/L)^3)*(O+t)*L,
#                   O = 0.9e-2 .. 0.11e-1
#                 ),
#                 Beeta = 0 .. b
#               );
#
# 'Simplification' isn't assigned to anything, and
# isn't subsequently used, so comment out
#
#  simplify
#  ( Mo1 = int
#          ( int
#            ( 2*shearmodulus*(1+Beeta*theeta/L-1/(1+Beeta*theeta/L)^3)*(O+t)*L,
#              O = 0.9e-2 .. 0.11e-1
#            ),
#            Beeta = 0 .. b
#          )
#  );
#
# Set up triple integral in 'Inert' form using
# Int() rather than int()
#
  Mo2 := 2*( Int
             ( Int
               ( Int
                 ( shearmodulus*((R+b+sin(Phi)*(A+tao))/R-1/((R+b+sin(Phi)*(A+tao))/R)^3)*((A+tao)^2*sin(Phi)+b*(A+tao))*L,
                   Phi = 0 .. (1/2)*Pi
                 ),
                 A = 0.9e-2 .. 0.11e-1
               ),
               tao = 0 .. t
             )
           );
#
# Evaluate the above numerically
#
  evalf(%);

  with(ListTools):
  L:=[a,b,c,d,e,f,g,h,i,j,k,l];
#
# Get the index of 'c' in L
#
  Search(c,L);
#
# Get the index of 'h' in L
# starting from the index of 'c'
# in L
  Search( h, L[Search(c,L)..-1]);
#
# Or possibly(?) return the list from
# index 1 to the first occurrence
# 'd' of
#
  L[1..Search(d, L)];
#
# Or possibly return the list from the
# first occurrence of 'c' to the first
# occurrence of 'g'
#
  L[Search(c, L)..Search(g, L)];

  restart:
  with(ExcelTools):

#
# generate an array of "points" and plot them
#
  points:= Array
           ( [ seq
               ( evalf~
                 ( [ cos(Pi*T/40),
                     sin(Pi*T/40),
                     T/40
                   ]
                 ),
                 T=0..40
               )
             ]
           ):
  plots:- pointplot3d
          ( points,
            connect=true,
            color=red
          );
#
# Save these points to an Excel file. NB
# OP will have to change the file path/name
# to something appropriate for his/her
# macin
#
  Export
  ( points,
    "C:/Users/TomLeslie/Desktop/testData.xlsx",
    "Sheet1",
    "B2"
  );

#
# Restart everything and load the
# ExcelTools() package
#
  restart;
  with(ExcelTools):
#
# Import the data from the Excel file
# and assign it to the variable 'points
#
  points:= Import
           ( "C:/Users/TomLeslie/Desktop/testData.xlsx",
             "Sheet1",
             "B2"
           ):
#
# Now plot the data in the variable
# 'points'
#
   plots:- pointplot3d
           ( points,
             connect=true,
             color=red
           );

  restart;
  Expr1:=(2*abs(x-2)+abs(x+1)-3)/(abs(x^2-4)+x-2);
#
# No problems with basic 'limit()' commands
#
  limit( Expr1, x=2, right);
  limit( Expr1, x=2, left);

#
# On the other hand the limitTutor  only
# finds one of these. The second reported
# limit is wrong!!!!
#
  Student[Calculus1]:-LimitTutor( Expr1, 2, 'right');
  Student[Calculus1]:-LimitTutor( Expr1, 2, 'left');

#
# Easier to apprecitate these limits by
# converting the original expression to
# 'piecewise'
#
  Expr2:=convert(Expr1, piecewise);
  plot(Expr2, x=-10..10, discont=true);
#
# Again no problems with basic 'limit()' commands
#
  limit( Expr2, x=2, right);
  limit( Expr2, x=2, left);

#
# On the other hand the limitTutor() now fails on
# both limits
#
  Student[Calculus1]:-LimitTutor( Expr2, 2, 'right');

  Student[Calculus1]:-LimitTutor( Expr2, 2, 'left');

  restart;
  a:= plot(1/x, x=0..1);
  b:=plottools:-getdata(a)[3][]:
#
# op(3,a) will return the 'view' setting
# which has (presumably) been determined
# by 'smartview'
#
  op(3,a);
#
# Hence can generate the data 'actually'
# being used by filtering the points in
# the getdata() matrix to leave only
# those within the view range
#
  M:= Matrix( [ seq
                ( `if`
                  ( `and`( b[j,1]>=op([3,1,1],a),
                           b[j,1]<=op([3,1,2],a),
                           b[j,2]>=op([3,2,1],a),
                           b[j,2]<=op([3,2,2],a)
                         ),
                     b[j, 1..2],
                     NULL
                  ),
                  j=1..237
                )
             ],
             scan=columns
           ):
#
# Determine ranges being plotted
#
  xRange:=min(M[..,1])..max(M[..,1]);
  yRange:=min(M[..,2])..max(M[..,2]);
#
# Check that the "revised" plot is the
# same as that plotted originally. Note
# that one has to adjust the 'view'
# option
#
  plot( M[..,1], M[..,2],
        view=[op([3,1],a), op([3,2],a)]
      );