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

There are 3 different kinds of plot3d-structures (see ?plot,structure), so the answer depends on your type of plot3d, as illustrated here.

#GRID
plot3d(x^2+y^2,x=-1..1,y=-1..1,numpoints=9);
p1:=%:
op([1,0],p1);
M1:=op(indets(p1,Array));
M1[1,1];
#MESH
plot3d(x^2+y^2,x=-1..1,y=-1..x,numpoints=9);
p2:=%:
op([1,0],p2);
M2:=op(indets(p2,Array));
M2[1,1];
#ISOSURFACE
plots:-implicitplot3d(z=x^2+y^2,x=-1..1,y=-1..1,z=0..2,grid=[5,5,5]);
p3:=%:
op([1,0],p3);
M3:=op(indets(p3,Array));
M3[1,1,1];

Preben Alsholm

You may try the following.

sys:={diff(cm(x,t),t)=Dm*diff(cm(x,t),x,x)-Vm*diff(cm(x,t),x)-pm*(cm(x,t)-cim(x,t)),diff(cim(x,t),t)=pim*(cm(x,t)-cim(x,t))}:
param1:={Vm=100,t0=0.0104167,a=5800,L=80};
conds:={cm(0,t)-Dm/Vm*D[1](cm)(0,t)=piecewise(t<t0,a,0),cm(L,t)+Dm/Dm*D[1](cm)(L,t)=0,cm(x,0)=0,cim(x,0)=0}:
PDE:=proc(p,pi,d) local param,res,RES,CM,CIM;
param:=param1 union {pm=p,pim=pi,Dm=d};
res:=pdsolve(eval(sys,param),eval(conds,param),numeric,time=t,range=0..80,timestep=1e-4);
RES:=res:-value([cm(x,t),cim(x,t)],output=listprocedure);
CM,CIM:=op(map2(subs,RES,[cm(x,t),cim(x,t)]));
proc(tt) local CMx; CMx:=fdiff(CM,[1],[40,tt]); eval(CM(40,tt)-Dm/Vm*CMx(tt),param) end proc;
end proc:

#Assuming that your (t,h)-points are in a list of lists called pkt:

F:=proc(p,pi,d) local H; H:=PDE(p,pi,d); add((H(pkt[k,1])-pkt[k,2])^2,k=1..80) end proc:

#Now using Sergey Moiseev's package DirectSearch (available from the Maple Application Center:

with(DirectSearch):
Search(F,initialpoint=[1,1,100],tolerances=.001);

The result is

[10.6213322849071, 5.70381751302765, 325.326647325447]

with an F-value of 

973.4

which is slightly better than the F-value of the results in the pdf-file you refer to:

984.7

I didn't have any success with the Optimization package in Maple.

Preben Alsholm

Another road to the final solution:

u:=sin(x)+sin(x+Pi/3):
expand(u);               
convert(%,phaseamp,x);

Preben Alsholm

An idea also used in a different thread started by Alec Mihailov:

S := Sum(ln(x)/(x^2+3*x+2),x=1..infinity):
Z :=k-> Sum(ln(x)/x^k,x=1..infinity):
value(Z(k));
                                 -Zeta(1, k)
asympt(ln(x)/(x^2+3*x+2),x);
                ln(x)   3 ln(x)   7 ln(x)   15 ln(x)    /1 \
                ----- - ------- + ------- - -------- + O|--|
                  2        3         4          5       | 6|
                 x        x         x          x        \x /
factor(combine(S-Z(2)+3*Z(3)));
                         infinity                  
                          -----                    
                           \                       
                            )     ln(x) (7 x + 6)  
                           /     ------------------
                          -----                   3
                          x = 1  (x + 2) (x + 1) x
evalf(%);
                                0.2586129327
%+evalf(value(Z(2)-3*Z(3)));
                                0.6017824583

Preben Alsholm