@acer Very helpful.  I see the difference between the "doable," and the "practicality." 

"...or they vary significantly in scale." 

I believe this is what's going on.  Some of those terms in the f functions are 10E-6 smaller than others, yet those terms contain the detail I want to calculate.  The Order (1) terms effectively cancel.  But I suspect the software has to keep track of terms like 1.00000010000 with at least that much accuracy.  As mentioned, I had been spending months on trying to get solutions from Engineering Equations Solver (EES, from fchart software), and I started seeing serious looking errors appearing.  The version of EES I used claims to carry calculations to 32 decimal places, so I don't think it's a round off problem, but something more basic, as you suggest.  Then again, maybe my old version of EES simply has a bug.  

I don't know what it means to treat it as an optimization problem.  Do you mean that in some way the demand for accuracy is tempered by a sufficient approach that is more practical?  Thanks for your suggestions. 

But thanks much for your suggestions.  

@ Venkit, I hope the improvements corrected those problems.  You might have a good suggestion to write a NR code myself.  

@Preben Alsholm 

Sorry I wasn't explicit enough.

I'm considering using fsolve for a simultaneous set of nonlinear algebraic equations containing hundreds of sine and cosine functions of the variables.

Interesting to me that the reference manuals do not specify any of the properties of their solver command (Nsolve, FindRoot, etc.), such as limitations on the nature of the arguments, maximum number of equations, maximum number of characters, etc. Am I to assume that since there are no stated limitations, there are none?

Important to me is whether I could define coefficients of the terms in the equations that are in turn expressed as a function of both independent variables and numerical constants. All the examples I see contain numeric coefficients.

For instance, here's what I'd like to do.

Independent Equations:

1) f1(X1, X2, X3, X4) = 0

2) C1*X1*Sin(X2) - C2*f2 = 0

3) C3*X1*(C4 – X4^2) – C4*f3 = 0

4) C5*X1*(Sin(X2) + Sin(X3)) – C5 = 0

where,

f2 = f2(X1, X2, X3, X4)

f3 = f3(X1, X2, X3, X4)

Each dependent function f is a sum of 22 components, with each component containing dozens of terms involving the four independent variables and other numerical parameters of the problem. Thus, up to a thousand or so total number of terms.

Can I simply define the f (dependent) functions separately from the fsolve command or do I need to convert the definitions of the f's into equations by moving everything to the left side and equating to zero?  That will give me a couple hundred equations.  Can Maple handle that?

Thanks for any comments and best regards.