restart;
#
# Matrix of dependent variables
#
  X:= Matrix( 2,1, [x__1(t), x__2(t)]);
#
# "Coefficient" matrix
#
  Mcoeff:= Matrix(2,2, [[1, -1], [12,-8]]);
#
# "Scalar addition" matrix
#
  b:=Matrix(2,1, [-3,4]);
#
# Define ode system in terms of above matrices
#
  ode:= diff(X,t) = Mcoeff.X + b;
#
# Solve ode system
#
  dsolve(ode);

  restart:
  interface(version);
  

  R:= 7.5;
  f:= z->z;
  a:= [1,2,3,4,5];
  b:= [2,1,1,2,4];
  g:= z-> CurveFitting:-Spline(a, b, z, degree=1);

#
# Example 1
#
  Int( Int( Int( arctan(f(r)),r=RR..R)/RR, RR), RR=0..R);
  evala(value(%));
#
# Example 2
#
  Int( Int( Int( g(r),r=RR..R)/RR, RR), RR=0..R);
  evala(value(%));
#
# Example 3
#
  Int( Int( Int( g(r)+arctan(f(r)),r=RR..R)/RR, RR), RR=0..R);
  evala(value(%));

  restart:
  interface(version);
  

  R:= 7.5;
  f:= z->z;
  a:= [1,2,3,4,5];
  b:= [2,1,1,2,4];
  g:= z-> CurveFitting:-Spline(a, b, z, degree=1);

#
# Example 1
#
  Int( Int( Int( arctan(f(r)),r=RR..R)/RR, RR), RR=0..R);
  evala(value(%));
#
# Example 2
#
  Int( Int( Int( g(r),r=RR..R)/RR, RR), RR=0..R);
  evala(value(%));
#
# Example 3
#
  Int( Int( Int( g(r)+arctan(f(r)),r=RR..R)/RR, RR), RR=0..R);
  evala(value(%));

Error, (in Handlers:-hasPiecewise) mismatched multiple assignment of 2 variables on the left side and 1 value on the right side

  restart;
  JointNum:= 3:
  alpha:= table( [seq(j=0, j=1..JointNum)]):
#
# test
#
  alpha[1];
  alpha[2];
  alpha[3];

  restart;
  JointNum := 3:
  seq( (alpha[j]:=0), j=1..JointNum):
#
# test
#
  alpha[1];
  alpha[2];
  alpha[3];

  restart;
  JointNum:= 3:
  seq( assign(alpha[j], 0), j=1..JointNum):
#
# test
#
  alpha[1];
  alpha[2];
  alpha[3];

  restart;
  JointNum:= 3:
  alpha:= Array(1..3, fill=0):
#
# test
#
  alpha[1];
  alpha[2];
  alpha[3];

  A:= u[1,2]+u[2]*u[1,1,1];
  LT:= (expr, term)-> select(has, expr, term);
  LT(A, u[1,1,1]);
  LT(A, u[2]);
  LT(A, u[1,2]);

  restart;
  gamma__1:= 2:
  k:=gamma__1*(100-x)/100-x*0.05*(gamma__1*(100-x)/100):
#
# Ensure that 'M' is only defined when k>0
#
  M:= 41539.42878*piecewise( k>0,
                             0.010+0.0001*((k+2*gamma__1)/2)^2,
                             undefined
                           ):
#
# Observation, k>0 when x<20 or x>100
#
  plot(k, x=-200..200);
  plot( abs(M), x=-200..200);

  restart;
  RA:= rtable(1 .. 5, 1 .. 5):

  ik:= 1:

  s:= proc()
           global ik, RA;
           RA[ik, 1]:= 1;
           ik:= ik + 1;
           RA[ik, 2]:= 2;
           ik:= ik + 1;
      end proc:

  t:= proc()
           global ik, RA;
           RA[ik, 3]:= 3;
           ik:= ik + 1;
           RA[ik, 4]:= 4;
           ik:= ik + 1;
       end proc:

  s():
  t():
  RA;
  ik;

  restart;
  expr:=u[1,1]*u[2,3]+u[1,1,2];
  (`[]`@ op)~(indets(expr, indexed));

  restart;
  with(geometry):
  with(plots):
  local O:
  _EnvHorizontalName := x: _EnvVerticalName := y:

  alias(coor = coordinates):
  ell := x^2/a^2+y^2/b^2 = 1:
  point(P1,a*cos(omega), b*sin(omega)):
  point(P2,a*cos(omega-(1/2)*Pi), b*sin(omega-(1/2)*Pi)):
  point(P3,a*cos(omega+(8/7)*Pi), b*sin(omega+(8/7)*Pi)):
  point(P4,a*cos(omega+5*Pi*(1/2)), b*sin(omega+5*Pi*(1/2))):
  a := 5: b := 3: omega := (1/5)*Pi:
  Ell := implicitplot(ell, x = -a .. a, y = -b .. b, color = red):
  dr := draw([seq(P || k, k = 1 .. 4)], axes = normal, printtext = true):

  for i from 1 to 4 do
      tgP||i := x*coor(P||i)[1]/a^2+y*coor(P||i)[2]/b^2 = 1
  od:
  poly := Matrix([coor(P1), coor(P2), coor(P3), coor(P4)]):
  Quadri := polygonplot(poly, axes = normal, color = "DarkGreen", transparency = .8):

  with(combinat): with(ListTools):
  L := [1, 2, 3, 4]:
#
# following for loop seems a bit redundant??
# Comment out for now
#
#  for i from 1 to 4 do
#      Rotate(L, i)[1]
#  od:

#
# Fix syntax in this loop
#
  for i to 4 do
      solve({tgP || (Rotate(L, i)[1]), tgP || i});
      point(S || i, subs(%, x), subs(%, y));
      coor(S || i)
  end do:

#
# Following code repeats what above 'for' loop
# does?? Why have it. Comment out for now
#

#  solve({tgP1, tgP2}, {x, y}):
#  point(S1, subs(%, x), subs(%, y)):
#  coor(S1):
#  solve({tgP2, tgP3}, {x, y}):
#  point(S2, subs(%, x), subs(%, y)):
#  coor(S2):
#  solve({tgP3, tgP4}, {x, y}):
#  point(S3, subs(%, x), subs(%, y)):
#  coor(S3):
#  solve({tgP1, tgP4}, {x, y}):
#  point(S4, subs(%, x), subs(%, y)):
#  coor(S4):
  
  poly:= Matrix([coor(S1), coor(S2), coor(S3), coor(S4)]):
  Quadri2:= polygonplot(poly, axes = normal, color = "DarkGreen", transparency = .9):
  #dr2:= draw(seq(S||k,k =1..4), axes = normal, printtext = true):
  line(diag13, [S1, S3]):
  line(diag24, [S2, S4]):
  midpoint(M1, S1, S3):
  midpoint(M2, S4, S2):
  line(Lm, [M1, M2]):

  dr2 := draw( [S1, S2, S3, S4, M1, M2, Lm(color = black), diag13, diag24],
               axes = normal,
               printtext = true
             ):

  for i from 1 to 4 do
      TgP||i := implicitplot( tgP||i,
                              x = -a-5 .. a+5,
                              y = -b-5 .. b+5,
                              color = blue
                            )
  od:

  display( [Ell, seq(TgP||i,i=1..4), Quadri, Quadri2, dr, dr2],
           view = [-a-5 .. a+3, -b-2 .. b+2],
           scaling = constrained,
           size = [700, 700]
         );
  
 

  restart;
  Chebyshev_approx:= proc(f::procedure, N::numeric)
                          local F:=0, a:= Array(0..N), n;
                          for n from 0 to N do
                              if   n=0
                              then a[n]:=1/Pi*evalf( Int(1/sqrt(1-x^2)*f(x), x=-1..1) );
                              else a[n]:=2/Pi*evalf( Int(1/sqrt(1-x^2)*f(x)*ChebyshevT(n , x), x=-1..1) );
                              fi;
                              F:=F+a[n]*ChebyshevT(n, x);
                          od;
                          return eval(a), F;
                    end proc:
  A, F:=Chebyshev_approx( x->sin(x), 6):
  A;
  F;

  restart:
  fName:="C:/Users/TomLeslie/Desktop/data.txt":

  myRead:= proc( myFile )
                 local j:=1,
                       l,
                       line;
                 while true do
                       line:=readline( myFile );
                       if   line <> 0
                       then l[j]:=[parse(line)];
                            j:=j+1;
                       else break;
                       fi;
                od:
                return Matrix([entries(l, 'nolist')])
         end proc:

   M:=myRead( fName)

   restart:

  b1:=1.41:d:=0.5/1:xi:=0.1:ea:=0.5:ra:=2:

  L:= unapply
      ( piecewise
        ( d<=z,
          1-2*xi*(cos((2*3.14)*((z-d)*(1/2))-1/4)-(7/100)*cos((32*3.14)*(z-d-1/2))),
          z<=d+1,
          1-2*xi*(cos((2*3.14)*((z-d)*(1/2))-1/4)-(7/100)*cos((32*3.14)*(z-d-1/2))),
          1
        ),
        z
      );
#
# OP seems interested in L(0.71) - not sure why
#
  L(0.71);
#
# OP also seems interested in the value of z fo
# which L(z)=1 (again not sure why)
#
  ans:=fsolve( L(z)=1,z);

#
# Define PDE and boundary/initial conditions
#
  PDE1 :=ra*(diff(f(x,t),t))=+b1*(1+ea*cos(t))+(1/(L(z)^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 for z=0.71 and obtain the value of 'f'
# at x=0.71, t=1.12
#
  sol:=pdsolve( eval( PDE1,
                      z = 0.71
                    ),
                IBC,
                numeric
              );
  sol:-value( f(x,t))(0.71, 1.12 );
#
# Idle curiosity plot f(x,t) as a function of x, t.
#
  sol:-plot3d( f(x,t), x=0..1, t=0..10);

#
# Solve the PDE for L(z)=1 and obtain the value of 'f'
# at x=0.71, t=1.12
#
  sol:=pdsolve( eval( PDE1,
                      z = fsolve( L(z)=1 )
                    ),
                IBC,
                numeric
              );
  sol:-value(f(x,t))(0.71, 1.12);
#
# Idle curiosity plot f(x,t) as a function of x, t.
#
  sol:-plot3d( f(x,t), x=0..1, t=0..10);