.

You have 44 equations in 41 variables and Maple would tell you, that in
your case it is inconsistent.

Denote M you matrix matcosdir and v = Fexternas+Rvinculo.

Then you use LinearAlgebra[LinearSolve] to solve the linear equations 
M . x = - v, your 'Normales' then is that x (if there is some, there
may be none and there may be lower-dimensional spaces as solution and
there may be exactly 1).

  with(LinearAlgebra);
  LinearSolve(M,-v);

    Error, (in LinearAlgebra:-LA_Main:-LinearSolve) inconsistent system

Actually you can check it (theoretical):

  interface(rtablesize=50);      # displays larger matrices
  ReducedRowEchelonForm(<M|-v>): # cf. help + Math books (help is vague)
  map(fnormal,%,3):              # you only have 3 Digits
  simplify(%,zero): identify(%);

or 

  GaussianElimination(M);

Your matrix then is zero in the last lines, while your v is not:
it is like solving 0*x=1 for x.


PS: your sheet is another example of the incredible uglyness of
the 'standard' worksheets, it does not even allow to save in a
reasonable format as *.wms ... what a shame for Maple.

This is why I usually refuse too look at *.mw (which unfortunately is
used by beginners).

The function should be neither too flat not too steep, since
you divide by f'(x[n]) and actually a code would check that.

Beyond that the theory says: do not start at inflection points.
But practical ...

So an answer would be: there is no guarantee that it works
and if not ... try another method ...

"Numerical Recipies" shows some, Google may give you
links for some downloads of that book.

PS: you should not use fdiff here, try D(f).

Well, I think there is no whole or circular argument anymore, but
it will depend on what you are using as defintion and properties
for exp(x) and/or a^x.

1. Version: like I did above, l'Hospital filling up arguments after
Alex's comments. The only thing missing there is exp(x) ' = exp(x).
That follows if you know the (absolute convergent) series for exp
and you know/accept differentiation of such series, since for term
n you have ( x^(n+1)/(n+1)! )' = (n+1)*x^n/(n+1)! = x^n/n!.

2. Version: if you are willing to see the task as difference version
for f'(x) in x=0, f:=a^x, you can even omit the dubious 'Hospital
and use the argumentation given there after that.

3. Version: I looked that up in the nice little book Otto Forster "Analysis1".

He introduces exp as series (and by quotient criterion it is seen absolute
convergent), then gives remainder estimate by geometric series, (10 lines)
error cutting after the N-th term <= 2*|x|^(N+1)/(N+1)!, for |x| <= 1+N/2.

Then |exp(x) - (1+x)| <= x^2 for |x |<= 1 + N/2, N=1 using that estimate.

Now dividing by 0 < |x| < 1 gives what you need

|  e^x - 1     |   | e^x - (1+x) |
|  ------- - 1 | = | ----------- |  <= |x|
|     x        |   |     x       |
 

  '[numer(f1),denom(f1)]': 'eval(%,h=0)': '%' =%;

                [numer(f1), denom(f1)]|      = [0, 0]
                                      |h = 0

  '[diff(numer(f1),h),diff(denom(f1),h)]': 'eval(%,h=0)': '%' =%;

           d             d            |
          [-- numer(f1), -- denom(f1)]|      = [ln(a), 1]
           dh            dh           |h = 0

so it is ln(a) by l'Hospital's rule.

Over time I simplify got used to the fact, that Maple does not always
allow to input things, as one would it expect from typing Mathematics.

Even if it mainly allows much intuitive handling - much better then
what I have seen from Mathematica (where I find it quite odd to write
Exp[x] for exp(x) (exp is lower case in Math and functions have round
brackets - besides the uncomfortable handling on European keyboards).
Or Pari or others (but my experience is rather limited on other CAS). 

And that handling for me was on major point to decide for Maple (after
using Mupad, which handled similar).

Part of 'exceptions' are indexed variables, so I simply try to avoid
them as far as possible.

In many cases x(n) instead of x[n] would do it as well (and for Math
it is the same, indexing more or less just means to use functions on
integer intervals)


Now for your problem (where the 'trick' with x(0) would not work either).

If you do not use the arrow notation but 'unapply' then it works and
you even can see, how Maple beats it into its needed shape:

  42;                   # write down the function expression
  y:=unapply(%, x[0]);  # make it a function by telling it the variable(s)

                            y := x_0 -> 42

One sees, that Maple constructs a function, but replaces x[0] by a symbol
(and thus silently solved its problem with x[0] internally).

Which mathematically is correct: by using the Math notation x[0] has to
be understood as bound variable (in some domain), but of course that does
not depend on the name of the variable.

Now you can use it, it takes indexed variables as input:

  'y( x[0] )': '%'=%;
                             y(x[0]) = 42




As others already mentioned neither x[0] nor x(0) are considered as symbols
(and that's what the error message for the intially entered command wants
to tell us). 

  x[0]; whattype(%); type(%%,symbol);

                               indexed
                                false


  x(0); whattype(%); type(%%,symbol);

                               function
                                false

However that's not what I would expect and my guessing is, that here the
type of x is given, not that of x[0] or x(0).

At least the help shows type( f(x,y,z), function ) as 'true', while in
Math I would say type( f, function ) = true.

But for all of the following ways one gets 'false':

 type( sin, function );

 X:= x -> x; type( X, function );

 x; X:= unapply(%,x); type( X, function );


However using 'type( %, procedure )' gives 'true'.


Which is not quite consistent with intuitive Math notation ... but
usually I can live with those handicaps ...

If computing diagonalizations for symmetric matrices: does Maple go through Eigenvalues/Eigenvectors only (needs roots of polynomials of may be higher order) or does it use other ways (have not tested it beyond dim=3 or 4 in concrete examples)?

try to either use 'formatted' in posting or just put it in a file (prefered: *.mws) and upload using the almsot invisible green arrow at the end of the menu bar of the board editor (visible, if logged in)

if you would use common notations, then your f = 1/(cos(x)^2),
so at Pi/2 ... hence you should not write down the integral

something like

for k  from 1 to 100 by  5 do
  cat(convert(k,string),".jpg");
en do;

Digits:=14;

numQ := proc(y, t1, t2, t3, t4)
local lnConst, e1, e2, p, lnc0, c0, c1, c2, aleph, q, 
  result, t12, t13, t15, t19, t22, t6, t7, t9, lnK5, ln_aleph;
  lnConst := -4681.92777205978276;
  e1 := 9.10323944425767683;
  e2 := 6.56420567208544207;
  p := 945.427297455447429;
  lnc0 := -5.34824719351398036;
  c1 := 16.5500650823099687;
  c2 := .763123640889490148e-1;
  t6 := .217003416872837899e-1*t3;
  t7 := .77256843031344e-1*t4;
  t9 := exp(-.22850024468343e-1*t1+t6-t7+e1);
  t12 := .908028439918035485*t3;
  t13 := .184630769011192971e-2*t4;
  t15 := exp(-.74140364919166*t2+t12+t13+e2);
  t19 := exp(.228500244683431746e-1*t1+t6-t7+e1);
  t22 := exp(.741403649191661292*t2+t12+t13+e2);
  aleph := y*t9+y*t15+t19+t22;
  ln_aleph := ln(aleph);
  lnK5 := lnc0+ln(y)*c1+ln_aleph*p-t4*t4-t2*t2-t3*t3-t1*t1-
    47.7572820906515031*ln_aleph*ln_aleph-ln(-y+1.0)*c2;
  q := lnConst+lnK5;
  result := exp(q);
  return result
end proc

# plot an example
plot('numQ'(0.1,tau,tau,tau,tau), tau= -2 .. 2, title="y=0.1");

The original task should be done by integrating Q, t in IR^4,
and y in ]0 ... 1[, for which t in ]-6 .. 6[^4 might be enough,
may be up to 8 is better. Have not checked it. And not even
tried to think about error estimates seriously. But claim it
is OK ... well, I am not forced to publish ...

  gc();
  infolevel[`evalf/int`]:=5:
  st:=time():
  
  yRange:= 0 .. 1;
  tRange:= -8 .. 8;
  relError:= 1.0 * 10^(-3);
  
  [t1 = -infinity..infinity, t2 = -infinity..infinity, 
    t3 = -infinity..infinity, t4 = -infinity..infinity]:
  subs(-infinity= op(1,tRange), infinity= op(-1,tRange),%):
  Z1:=op(%):
  
  Int('numQ6', [yRange, tRange$4], method=_cuhre, epsilon = relError);
  evalf(%);
  
  `seconds needed` = time()-st;
  infolevel[`evalf/int`]:=0:

It needs ~ 13 sec on my PC and ~ 10 min for relError:= 1.0 * 10^(-6),
where in the latter case 57704829 evaluation points have been used,
6*10^7 function evaluations ... which is a lot.

No better timings if 'compile' is used, the integrator does not like it.


Just a summary on a fairly non-trivial task:

The above comes through processing the original task symbolically
for z and the last 4 variables, hence only 5 variables of the 10
survive. After that a 4-dim transform is found, which makes it an
integral ~ f(y,t)*pdfN(4,t), the later a 4-dim normal pdf with all
correlations = 0, hence centered at 0 (otherwise ~ (-10,100,0,0)
should be choosen, depending on choosing integration order) and it
dominates f for large t's. And finally all can be brought into the
posted procedure using codegen/optimze after resolving overflow issues.


That is my last post on that topic & thread.

Cheers

To me it looks as a system of quadrics and generically they intersect
discrete due to dim = #equations. But the only thing I can 'see' or
guess is that the first 5 are equivalent to "Mittelpunktsflächen"
(= "centered hypersurfaces" ?) and the others seem parabolic. Not sure,
but if one would have the normal form for them, then it may be clearer.

The reason for the question might certainly help ... sqrt(3) smells
a bit like something based on Algebra or Number Theory.

And if the solution is not a discrete set, then it might matter over
which ring this is considered (not that I want to solve it ... never
used that corner of Maple).

So what are your findings so far?

Do you get 'reasonable' results, if you use ~ 0.0001398 as result? I think, it could/should be of that size.

showstat(`evalf/exp`) or showstat(`exp`)shows the code for floats or symbolics