of this question is GI. After all, if you provide GI,  what do you expect to happen - a "sensible" answer from a garbage question?.

So let's unpick what is happening when you issue the command

diff(x-1=0,x)

the diff() command will cause both sides of the equation to be differentiated wrt to the variable 'x'. In other words you are making the statement

 diff(x,x)+diff(-1,x)=diff(0,x)

which of course will result in

1-0=0

Just out of curiosity, what answer did you expect? Because for this example of GI, I honestly have no idea!

If you really want to go achieve an award for GI, then forget all the complicated stuff like differentiation and just try typing

123456=0

in Maple. Now exactly what do you think should happen when you do this?

I have made a number of changes in your worksheet. You should read the comments in the attached very carefully just to confirm that these changes can only "improve" the situation.

odeProb.mw

However the ODE system still has a left end-point singularity, and applying a mid-point method does help (much), which could still be a result of the singularity at eta=0. You should give very careful thought to why such a singularity might occur. If your equations represent a "real-world" system - ask yourself the question: Should there be a singularity at the left end point??

A (fairly common) way to produce such an effect, is simply that the equations have been incorrectly entered!

If you are convinced that the equations are (now) OK, then report back and we can consider possible strategies

is probably to use the display() command with the 'insequence=true' option as in the following code

  restart:
  with(plots):
  with(plottools):
  display( [ seq
             ( cone( [0, 0, 1+j], j, -1-j ),
               j=0.1..1, 0.05
             )
           ],
           style=surface,
           insequence=true
         );

to do this is shown below

  restart;
#
# Return the inert sum
#
  H1a := Sum(GAMMA(2*b - 1 + 2*n)*GAMMA(2*n - s)*(b - 1/2 + 2*n)/(GAMMA(2*b + 2*n + s)*GAMMA(2*n + 1)), n = 0 .. infinity);
#
# use the value() command to evaluate the sum (if possible)
#
  value(H1a);

the attached.

constrMat.mw

You should be aware that I can't verify this in Maple 15: but I don't anticipate any problems in execution. The only thing which definitely won't work is the matrix display control  command

interface(rtablesize=[M1*M2,10])

is to use the "inert" version of the cos() function, whihc would be written %cos()

Note that this representation will stay "inert" throughout the rest of your worksheet untila you decide to make it non-inert, by using the value() command. Consider the output of the code

  restart;
#
# Make the cos() function inert
#
  eqaux_2 := l__bb = L__aa0 + L__aa2 *%cos(2*theta - 4/3*Pi);
#
# Now make the cos() function "active" again
  value(%);

(I think!)

If I fix various typos, incorrect use of indexing etc stc, then I come up with the code

  restart;
  N:= 1:
  u(r,z):= add( lambda[1]^j*u[j](r,z),j=0..N);
  h:= (1+xi*z)*(1-eta[1]*(z-sigma)-(z-sigma)^n);
  Eq1:= (1/r)*diff((r*(diff(u(r, z),r))+lambda[1]*(2*diff(u(r, z),r)*diff(u(r, z),z))),r):
#
# NB sol[-1] does nothing - It is just a harmless
# initialisation which allows the current iteration
# of the loop code to rerieve an answer from the
# previous iteration (even for j=0, when the "previous"
# iteration does not exist)
#
  sol[-1]:=u[-1](r,z)=0:

  for j from 0 to N do
      cond:= D[1](u[j])(0,z) = 0,
             u[j](h,z)= 0;
      sol[j]:= pdsolve
               ( { eval
                   ( coeff(Eq1, lambda[1], j) = 0,
                     sol[j-1]
                   ),
                   cond
                 },
                 u[j](r,z)
               );
  od;

This will provide an explicit answer for u₀(r,z), but for u₁(r,z), pdsolve() returns (a fairly lengthy) expression containing integrals. I'm not convinced that the OP's "iterative solution" approach is going to be productive. See the attached

pdeProb.mw

It seems as if what you actually want is a "kronecker delta" function, where

kron(a,b)=1 if a=b, 0 otherwise

which is pretty trivial to write, as in

kron:=(n1, n2)-> `if`(n1=n2, 1, 0)

This would be OK if the passed arguments could be guaranteed to be integers, but in your case  it seems that at least one of these will be a "float", so it may be safer to define the same thing with a "tolerance" parameter, something like

tol:=1e-09:
kron:=(a1, a2)->`if`( abs(a1-a2)<tol, 1 ,0)

You are still likely to run into problems when plotting any function which contains calls to kron(). Maple's default plotting algorithm does not allow you to specify exactly the points at which the plot will be evaluated, so may "miss" the points at which the kron(0 functin is non-zero. It is possible to circumvent this by constructing your own sequence of point to be plotted. You might want toi consider the "toy" example in the code below

kron:=(n1, n2)->`if`(n1=n2,1,0);
plot( [seq( [j, j^2*(kron(j,1)+kron(j,2))], j=0..3, 0.01)])

Hope this helps

For the parameters of your problem, the evaluation of Ip produces a complex number. If I check with

Ip :=evalf~( (1/3)*(eval(C, [x = X-(1/3)*L0, L = L0])));
printf("%a", Ip);

this returns

[.5000000000, 1.118033988*I]

Not that the y-coordinate is imaginary

1. Installing the Maple add-in to your Office version allows you to call Maple functionality from within Excel. (Generally speaking this is a bit "clunky" but it works). It has nothing to do with accessing Excel data from within Maple
2. If you are having trouble accessing data in an Excel file from within Maple, there are a couple of possibilities
   1. Incorrect file path - try using a full path name to the Excel file, ie (on Windows) ExcelTools:-Import( "C:/path/path/path/Excelfile", "sheetname", "cellnames")
   2. Verify permissions - who has permission to read your Excel file? Are you certain that the "Maple-user" has read privileges to the Excel file?

As a test, you can try the following code

#
# Construct a random matrix
# 
  M:=LinearAlgebra:-RandomMatrix(8,8);
#
# Export this matrix to an Excel file. For
# test pruposes, I am using my default "desktop"
# on Windows. You wil have to change this to
# something appropriate for your machine. Also
# for test purposes I am sprecifying an Excel
# sheet name of "mySheet", and that M[1,1]
# should be in the Excel cell B2
#
  ExcelTools:-Export( M,
                      "C:/Users/TomLeslie/Desktop/test.xls",
                      "mySheet",
                      "B2"
                    ) 
#
# Restart everything
#
  restart;
  M:=ExcelTools:-Import( "C:/Users/TomLeslie/Desktop/test.xls",
                         "mySheet",
                         "B2:C3"
                       );

which writes a matrix to an Excel file, restats, then reads a portion of that matrix. This works perfectly.

In order to make this code run on you machine, you will have to change the filepath "C:/Users/TomLeslie/Desktop/" to somethng appropriate for yur machine

I don't think that the help is "misleading", but it is worth reading very carefully. The most important point to grasp is (snip from ?PositionVector, emphasis added),

The position Vector is a Cartesian Vector rooted at the origin, and has no mathematical meaning in non-Cartesian coordinates.

Thus when you enter the command

      v1:=VectorCalculus:-PositionVector([r,0,0], spherical);

[r,0,0] are the components in spherical polars, and these components will be "converted" to the equivalent cartesian components when the position vector is created. The positionVector() command always returns a vector in cartesian coordinates. For this specific example, since the polar angle is 0, the value of the azimuthal angle is irrelevant, and the resulting position vector will lie along the positive z-axis. This is confirmed with the command

VectorCalculus:-About(v1);

which returns

                     Matrix(4, 2, [ [`Type: `, `Position Vector`],
                                          [`Components: `, [0, 0, r]],
                                          [`Coordinates: `, cartesian],
                                          [`Root Point: `, [0, 0, 0]]
                                       ]
                              )

This behaviour becomes even more obvious if you look at a "general" point in spherical polars, with

v2:=VectorCalculus:-PositionVector([r,theta,phi], spherical);
VectorCalculus:-About(v2);

where the About() command returns

Matrix(4, 2, [[`Type: `, `Position Vector`],
                    [`Components: `, [r*sin(theta)*cos(phi), r*sin(theta)*sin(phi), r*cos(theta)]],
                    [`Coordinates: `, cartesian],
                    [`Root Point: `, [0, 0, 0]]
                 ]
           )

Not really clear wnat yur problem is but the following code will initialise (and extend) your Arrays more efficiently as weel as producing the DataFrames yu seem to want

  restart;
  with(ArrayTools):
  nOben:=10:
  nUnten:=3:
#
# Initialise two vectors
#
  InitialisierungDF:=Vector[column](nOben-nUnten+1, fill=oE):
  InitialisierungSpalte:=Vector[row](nOben-nUnten+1, i->n=i):
#
# The following will produce two identical Data Frames
#
# Why two?
#
  DF1:= DataFrame( Concatenate( 2, InitialisierungDF $ 4), 
                   columns = [ GK, PZP, PZM, PY],
                   rows = InitialisierungSpalte);  
  DF2:= DataFrame( Concatenate( 2, InitialisierungDF $ 4),
                   columns = [ GK, PZP, PZM, PY],
                   rows = InitialisierungSpalte);
  DF3:= DataFrame( Concatenate( 2, Concatenate( 1, InitialisierungDF, <oE>) $ 4),
                   columns = [ GK, PZP, PZM, PY],
                   rows = Concatenate(2, InitialisierungSpalte, RMSE)
                 );

Your problem is that the intgral is being evaluated before the change of variables can be performed. You can avoid this by making the integral "inert", using

IntegrationTools:-Change(Int(3*x*sqrt(x + 8), x), x = u^2 - 8);
value(%);
subs(u^2 = x + 8, %);

See

intInt.mw

which is a useful addition to Maple literature.

All of the code in this book is written with the assumption that the reader will be using 1-D input in worksheet mode, which is as WYSIWYG as you can get

Unfortunately, the current default for new Maple users is to use 2-D input in document mode. Entering (intended) 1-D code in a worksheet set up for 2-D input is guaranteed to cause chaos. As a smple illustration, consider typing

f(x):=x^2;

In a worksheet set up for 1-D input this is a simple assignment of the form name:=expression.  In particular the left-hand side of the above assignment is a simple 'name', which just happens to contain  the characters '(' and ')'

In a worksheet set up for 2-D input, the parser will interpret this as defining a function f() which accepts an argument 'x' and returns x^2. It is not a simple assignment, it is equivalent to writing f:=x->x^2.

Unfortunately for the OP, Articolo's book uses many simple assignments where the 'name' on the left-hand side contains parentheses

I have just painstakingly entered Articolo's code for the example 8.4.1 using 1-D input, and it runs perfectly, producing the graphs and animation of the book

artic.mw

  restart;
  with(plots):
  a:=1: k:=1/20: h(x,t):=0:
  A(t):=10: B(t):=20: kappa[1]:=1: kappa[2]:=0: kappa[3]:=1: kappa[4]:=0:
  b(t):=(kappa[3]*a*A(t)+A(t)*kappa[4]-B(t)*kappa[2])/(kappa[1]*kappa[3]*a+kappa[1]*kappa[4]-kappa[2]*kappa[3]);
  m(t):=(kappa[1]*B(t)-A(t)*kappa[3])/(kappa[1]*kappa[3]*a+kappa[1]*kappa[4]-kappa[2]*kappa[3]);
  s(x,t):=m(t)*x+b(t):
  s(x,0):=eval(subs(t=0, s(x,t)));
  lambda[n]:=(n*Pi/a)^2;
  X[n](x):= sqrt(2/a)*sin(n*Pi/a*x); X[m](x):=subs(n=m, X[n](x));
  w(x):=1:
  Int(X[n](x)*X[m](x)*w(x), x=0..a)=delta(n,m);
  diff(T[n](t),t)+k*lambda[n]*T[n](t)=Q[n](t);
  q(x,t):=h(x,t)-diff(s(x,t),t);
  Q[n](t):=Int(q(x,t)*X[n](t),x=0..a); Q[n](t):=expand(value(%)):
  Q[n](t):=simplify(factor(subs({sin(n*Pi)=0,cos(n*Pi)=(-1)^n}, Q[n](t))));
  Q[n](tau):=subs(t=tau, %):
  T1(t):=exp(-k*lambda[n]*t);
  G1(t,tau):=simplify(T1(t)/(subs(t=tau,T1(t))));
  T[n](t):=eval(C(n)*T1(t)+Int(G1(t,tau)*Q[n](tau),tau=0..t)); T[n](t):=value(%):
  v[n](x,t):=simplify(eval(T[n](t)*X[n](x)));
  v[n](x,t):=Sum(v[n](x,t),n=1..infinity);
  f(x):=60*x-50*x^2+10;
  C(n):=Int((f(x)-s(x,0))*X[n](x),x=0..a);C(n):=expand(value(%)):
  C(n):=simplify(subs({sin(n*Pi)=0, cos(n*Pi)=(-1)^n},C(n)));
  T[n](t):=(C(n)*T1(t)+int(G1(t,tau)*Q[n](tau),tau=0..t));
  v[n](x,t):=simplify(eval(T[n](t)*X[n](x)));
  u(x,t):=eval(Sum(v[n](x,t),n=1..infinity)+s(x,t));
  u(x,t):=s(x,t)+sum(v[n](x,t), n=1..3):
  animate( u(x,t), x=0..a, t=0..5, thickness=3);
  u(x,0):=subs(t=0,u(x,t)):u(x,1):=subs(t=1,u(x,t)):
  u(x,2):=subs(t=2,u(x,t)):u(x,3):=subs(t=3,u(x,t)):
  u(x,4):=subs(t=4,u(x,t)):u(x,5):=subs(t=5,u(x,t)):
  plot({u(x,0), u(x,1), u(x,2), u(x,3), u(x,4), u(x,5)}, x=0..a, thickness=5);

the problem is a lot simpler.

As Kitonumm has pointed out, for general (ie non-regular) tetrahedra, there may not even be an in-sphere.

However, a regular tetrahedron can be defined in terms of its "centre" and the in-sphere radius. Thhis makes generating the tetradrons and the in-sphere particularly simple. See the code below

  restart:
  with(geom3d):
  _EnvXName := x:
  _EnvYName := y:
  _EnvZName := z:
#
# Define the "centre" of the tetrahedron
# Note that [0,0,0] is just a kind of obvious
# choice - feel free to pick any other value
#
  point(o,0,0,0):
#
# Define the "in-radius" of the tetrahedron
#
  inRad:=2:
#
# Define the tetrahedron with the centre and
# in-radius specified above
#
  tetrahedron( tet, o, in_radius=inRad):
#
# Define the sphere with the centre and radius
# specified above
#
  sphere(sph, [o, inRad]):
#
# Draw the tetrahedron and in-sphere
#
  draw( [ tet(color=yellow, transparency=0.5),
          sph(color=green,  transparency=0.5)
        ]
      );
#
# Check the details of the tetrahedron and the
# in-sphere
#
  detail(tet);
  detail(sph);

which will produce the plot below (along with  other details about the tetrahedron and the in-sphere)