To see what is going on, modify the example in the help for Statistics[PolynomialFit]:

with(Statistics):
X := Vector([1, 2, 3, 4, 5, 6], datatype=float):
Y := Vector([2, 3, 4, 3.5, 5.8, 7], datatype=float):
PolynomialFit(5, X, Y,x);
Warning, there are zero degrees of freedom
                                                                   2
22.3000000000013934 - 45.8000000000009352 x + 36.2083333333324831 x

                          3                        4                         5
   - 12.5833333333326678 x  + 1.99166666666652214 x  - 0.116666666666656538 x
PolynomialFit(6, X, Y,x);
Error, (in Statistics:-PolynomialFit) the number of independent parameters is greater than the number of observations

When you ask for the polynomial to have degree 5 there are 6 coefficients to be determined, and you have 6 equations, one for each of the (x, y)-points. Thus most likely the coefficients are uniquely determined, and the resulting polynomial is an interpolating polynomial passing through all of the 6 points.

You could try

?PDEtools[declare]

As I understand the question, the maximum has to be found for each of 100 values of G.

Using GlobalSearch from Sergey Moiseev's package DirectSearch (available from Maple Application Center) with G = 0.15:

GlobalSearch(eval(TF,G=0.15), {B >= 0, 0 <= a, B <= 1, a <= 1}, maximize);

#Output a = 0, B = 0.7841 agrees with the plot:

                  
plot3d(eval(TF,G=0.15),a=0..1,B=0..1,axes=boxed);

The basic reason for the failure is that f(z) = exp(1/z) has an essential singularity at 0.

However, with
F:=sinh(sqrt(phi/nu+m^2)*h);
G:=sinh(sqrt(phi)*h/sqrt(nu));

you can use your result to conclude that F/G -> 1 as nu->0, so that the main contribution to F for small nu is G.

To improve further you could use that for large x  sinh(x) = exp(x)/2 and then

restart;
G:=(1/2)*exp(sqrt(phi)*h/sqrt(nu)):
F:=1/2*exp(sqrt(phi/nu+m^2)*h):
combine(F/G):
series(%,nu=0,2):
G*%;

Of course the number 2 in series can be changed.           

If you increase the value of Digits to 20 before your calculations:

Digits:=20:

then you get values that are close.  

The equations eqn2 and eqn3 of course mean that N1 and N3 are constants.

If you multiply eqn1 by 2*diff(theta(s),s) and then integrate w.r.t. s you get

diff(theta(s),s)^2 = (2/15)*(N3*sin(theta(s))+N1*cos(theta(s))) + C

where C is a constant of integration.

To solve this equation you would need to take the square root, thus you get two branches and you will need to switch between these.

The possibility of constant solution parts in the 1. order equation is not in the original problem. It was introduced by the multiplication by diff(theta(s),s).

The error message from Maple is not very telling.

fsolve is used for numerical solution, so the command would be

fsolve(eq2,phiS);

I tried

?undocumented

and this gave the help page 'Undocumented Protected Names'. Some of those are procedures.

The last two arguments to Do should be arguments to display, like this:

Do( %Plot0 = plots[display]( plottools[ellipse]([xpos,ypos],3,4),view=[-10..110,-10..110],axes=boxed ) );

Except for one extra parenthesis (which just created a syntax error) Maple 12, 13, and 14 computed the integral without problem:

restart;
with(VectorCalculus):
X := r -> sqrt(2/Pi)*alpha*exp(-alpha^2*(r.r)):
a_0:=a* <1/2*sqrt(3) , 1/2>:
r:=<x,y>:
assume(alpha>0);
about(alpha);
eq1:=int( X(r-a_0) , x=-infinity .. infinity);

The last assumption on a variable overwrites one made previously. If you want to make additional assumptions on a variable then use 'additionally'.

However, if you assume that alpha is positive, then you are assuming that it is in the real range from 0 to infinity (see the output from 'about(alpha)' ), so there is no need for an additional assumption about alpha being real.

Without the VectorCalculus package you need to be aware that the inner product by default works in the complex domain:

restart;
X := r -> sqrt(2/Pi)*alpha*exp(-alpha^2*(r.r)):
a_0:=a* <1/2*sqrt(3) , 1/2>:
r:=<x,y>:
r.r;
eq1:=int( X(r-a_0) , x=-infinity .. infinity):
simplify(eq1) assuming real;
#The following syntax can be used in Maple 14:
simplify(eq1) assuming real,alpha>0;

Preben Alsholm

evalm and &* was used with the matrix-structure ('m', not 'M'). With A and B defined to be Matrices (capital M) you can just do

M:=A.B;

Preben Alsholm

I guess you could use matrixplot from the plots package, but why should it be faster to create the matrices from f1, f2, etc. and then plot their sum in some way?

Preben Alsholm

The range is split incorrectly in the third version, so no wonder that you get a wrong result.

Simpson needs more intervals, so try:

restart;
E_field := conjugate((-.7147663039-.8729430992*I)*exp((9.123185068*I)*cos(phi))*cos(phi)+(0.3187576278e-1+0.3371906130e-1*I)*sin(phi)*exp((4.561592534*I)*(-1.*sin(phi)+cos(phi)))+(0.7062621752e-1+.1302618973*I)*cos(phi)*exp((2.000000000*10^(-11)*I)*(53.*cos(phi)+1.570796328*10^11))+(-0.3187599160e-1-0.3371906130e-1*I)*sin(phi)*exp((4.561592534*I)*(sin(phi)+cos(phi))))*sin(phi):
Student:-Calculus1:-ApproximateInt(E_field, phi = 0..2*Pi,method=simpson,partition=20);
evalf(Int(E_field, phi = 0..2*Pi));
evalf(Int(E_field, phi = 0..5000/5001/2*Pi) + Int(E_field, phi = 5000/5001/2*Pi..5000/5001*Pi) + Int(E_field, phi = 5000/5001*Pi..5000/5001*3/2*Pi) + Int(E_field, phi = 5000/5001*3/2*Pi..2*Pi));
evalf(Int(E_field, phi = 0..2*Pi, method = _d01akc));
#

#Or with more digits:
Student:-Calculus1:-ApproximateInt(E_field, phi = 0..2*Pi,method=simpson,partition=100);
evalf[20](Int(E_field, phi = 0..2*Pi));
evalf[20](Int(E_field, phi = 0..5000/5001/2*Pi) + Int(E_field, phi = 5000/5001/2*Pi..5000/5001*Pi) + Int(E_field, phi = 5000/5001*Pi..5000/5001*3/2*Pi) + Int(E_field, phi = 5000/5001*3/2*Pi..2*Pi));
evalf[15](Int(E_field, phi = 0..2*Pi, method = _d01akc,epsilon=1e-11));

Preben Alsholm

Doesn't modifying the sum Tr do the job?

I'm assuming that by Tr:=sum(f, i = k .. n);
you meant
Tr:=sum((yei-f(x1i,x2i,x3i,x4i))^2,i=k..n)
Modified:
Tr:=sum((yei-f(x1i,x2i,x3i,x4i))^2+(zei-f(x1i,x2i,x3i,x4i))^2,i=k..n);

Preben Alsholm

If you use pecewise to begin with there is no problem:

f:=(x,lambda)->piecewise(x*lambda<0,x*lambda,0);
plot3d(f(x,lambda),x=-2..2,lambda=-2..2, axes=boxed);
D[2](f);
f2:=diff(f(x,lambda),lambda);
plot3d(f2,x=-2..2,lambda=-2..2, axes=boxed);

Preben Alsholm