> N:=4:

> U:=Vector[column](N,symbol=u):

> A:=Matrix(N,N,(i,j)->-1):

> f:=exp(U^%T.A.U):

> Int(f,seq(U[i]=1..infinity,i=1..N));

            /infinity   /infinity   /infinity   / infinity                              
            |           |           |           |                                        
            |           |           |           |           
            |           |           |           |                                        
           /1          /1          / 1         / 1                                       

              exp(u[1] (-u[1] - u[2] - u[3] - u[4]) + u[2] (-u[1] - u[2] - u[3] - u[4])

                  + u[3] (-u[1] - u[2] - u[3] - u[4]) + u[4] (-u[1] - u[2] - u[3] - u[4]))

              du[1] du[2] du[3] du[4]

> value(%);

                            35   (1/2)   17            35   (1/2)       
                          - -- Pi      + -- exp(-16) + -- Pi      erf(4)
                            6            12            6                

By transforming by logarithm, I got 4 roots even with default Digits=10. I used Digits=15 below just to get a more useful result, and I used Digits=100 only for checking the results.

restart:
eq:=tan(4*(1.1575*10^12/(1-.5)^1.2-x^2)^(1/2)*10^(-6)*(1/2))
    /tan(4*(5.555*10^11/(1-.5)^1.2-x^2)^(1/2)*10^(-6)*(1/2))
  = -4*x^2*(5.555*10^11/(1-.5)^1.2-x^2)^(1/2)*(1.1575*10^12
    /(1-.5)^1.2-x^2)^(1/2)/(1.1575*10^12/(1-.5)^1.2-2*x^2)^2:

Digits:=15:
ans:=Student:-Calculus1:-Roots(eq,x = 1*10^5 .. 9*10^6);

            [                   5                     6                     6]
            [2.38988443815009 10 , 1.12925289978792 10 , 1.74437318531275 10 ]

neweq:=subs(x=10^(y),eq):
Digits:=15:
newans:=Student:-Calculus1:-Roots(neweq, y = log[10](1*10^5) .. log[10](9*10^6));

                 [5.37837690133009, 6.05279121445174, 6.21237849725823, 

                  6.24163940203245]

Digits:=100:
newans:=map(t->10^t,newans):
Digits:=15:
newans:=evalf(newans);


               [                   5                     6                     6  
               [2.38988443815006 10 , 1.12925289978792 10 , 1.63071661915857 10 , 

                                   6]
                1.74437318531274 10 ]

Digits:=100:
check:=map(t->eval((rhs-lhs)(eq),x=t),newans):
Digits:=15:
evalf(check);

           [                   -15                      -11                     -8    
           [2.23350561239722 10   , -1.15959577074243 10   , 5.57289204905341 10   I, 

                                  -15]
              -8.86871138407261 10   ]

Try

     plot(<CosAa|CosBb>)

I didn't try to Compile the last version (I'm not sure if that would handle the nested add in the return line, or have other problems).

Are these results accurate?

restart:
NN:=7054:
M := Matrix(NN, 60):
for i to 10 do M[1, i] := 1 end do:
y:=[0.8045, 0.1834, 0.0006]:
Mat_proc := proc(x)
  local i, j;
  global M;
  for i from 2 to 7054 do
    M[i, 1] := M[i-1, 2];
    M[i, 60] := M[i-1, 59];
    for j from 2 to 59 do
      M[i, j] := M[i-1, j] + x * (M[i-1, j-1] - 2*M[i-1, j] + M[i-1, j+1])
    end do
  end do;
return [seq((1/10)*add(M[7054, i+10*j], i = 1 .. 10), j = 1 .. 3)]
end proc:
CodeTools:-Usage(add((Mat_proc(0.3)[i]-y[i])^2, i = 1 .. 3));
memory used=374.11MiB, alloc change=218.58MiB, cpu time=4.17s, real time=4.18s
                                0.4419552140

Mat_proc_new := proc(x,M,NN)
  local i, j;
  for i from 2 to NN do
    M[i, 1] := M[i-1, 2];
    M[i, 60] := M[i-1, 59];
    for j from 2 to 59 do
      M[i, j] := M[i-1, j] + x * (M[i-1, j-1] - 2*M[i-1, j] + M[i-1, j+1])
    end do
  end do;
return add(((1/10)*add(M[NN, i+10*j], i = 1 .. 10)-y[j])^2, j = 1 .. 3);
end proc:
CodeTools:-Usage(Mat_proc_new(0.3,M,NN));
memory used=82.37MiB, alloc change=0 bytes, cpu time=1.40s, real time=1.41s
                                0.4419552140

with(Optimization):
CodeTools:-Usage( Minimize('Mat_proc_new'(x,M,NN), x=1e-3 .. 0.2) );
memory used=0.95GiB, alloc change=186.84MiB, cpu time=18.56s, real time=18.65s
                [0.247283044276117, [x = 0.0120636643459173]]

Mat_proc_newer := proc(x,M,NN,y)
  local i::integer, j::integer;
  for i from 2 to NN do
    M[i, 1] := M[i-1, 2];
    M[i, 60] := M[i-1, 59];
    for j from 2 to 59 do
      M[i, j] := M[i-1, j] + x * (M[i-1, j-1] - 2*M[i-1, j] + M[i-1, j+1])
    end do
  end do;
return add(((1/10)*add(M[NN, i+10*j], i = 1 .. 10)-y[j])^2, j = 1 .. 3);
end proc:
Mf:=Matrix(M,datatype=float[8]):
yf:=Vector(y,datatype=float[8]):
CodeTools:-Usage(Mat_proc_newer(0.3,Mf,NN,yf));
memory used=79.00MiB, alloc change=0 bytes, cpu time=1.23s, real time=1.24s
                              0.441955214165890

CodeTools:-Usage( Minimize(x->Mat_proc_newer(x,Mf,NN,yf), 1e-3 .. 0.2) );
memory used=76.78MiB, alloc change=0 bytes, cpu time=4.40s, real time=4.46s
                  [0.247283044276117, [0.0120636643459173]]

#plot('Mat_proc_newer'(x,Mf,NN,yf), x=1e-3 .. 0.5, view=0..0.5);

Why not simply use the LinearAlgebra:-LinearSolve command?

Anyway, you can do substitutions into columns of a Matrix using usual rtable indexing and assignment.

G:=Matrix(2,2,[[a,b],[c,d]]);

                                      [a  b]
                                 G := [    ]
                                      [c  d]

V:=Vector(2,[e,f]);

                                       [e]
                                  V := [ ]
                                       [f]

cG:=copy(G):

cG[..,2]:=V:

cG; # after substitution

                                   [a  e]
                                   [    ]
                                   [c  f]

G; # unchanged

                                   [a  b]
                                   [    ]
                                   [c  d]

How about this?

restart:

myf := Int(exp(-(x-mu)^2/2/sigma^2)/sigma/sqrt(2*Pi)/(1+exp(-x)),
           x=-infinity..0, method=_d01amc)
     + Int(exp(-(x-mu)^2/2/sigma^2)/sigma/sqrt(2*Pi)/(1+exp(-x)),
           x=0..infinity, method=_d01amc):

plot3d(myf, mu=-5..5, sigma=0.5..10, axes=box);

It depends upon the nature of the problem (ie. univariate vs multivariate) and the setting of Digits.

evalf(Int(...)) has quite a lot of userinfo messages, at various levels of verbosity. Try the following examples as an illustration of some univariate cases. Notice how the print-out mostly seems to corroborate the claims about default methods made on the help page, in the section, "Outline of the Numerical Integration Polyalgorithm (1-D Integrals)."

I see a difference, though, running this in Maple 15.01. The help page states that the switchover between NAG methods and _CCquad happens when Digits becomes greater than evalhf(Digits), ie. when Digits=16 or higher. But the evidence below seems to indicate that it happens at a lower cutoff. I wonder whether some guard digits are added to Digits internally, before the decision is made. Or perhaps the decision process is more complicated.

restart:
infolevel[`evalf/int`]:=1:

Digits:=10:
evalf(Int(exp(-x)/x^2,x=1..100)); # first tries _d01ajc from NAG
evalf(Int(exp(-x)/x^2,x=1..infinity)); # first tries _d01amc from NAG

Digits:=12:
evalf(Int(exp(-x)/x^2,x=1..100)); # first tries _CCquad, moves on to _Dexp
evalf(Int(exp(-x)/x^2,x=1..infinity)); # first tries _CCquad, moves on to _Dexp

For your question 2)

Download recur.mw

There is no such thing as Maple 15.1, but you should be able to get the update to 15.01 by following this link.

`zip` and elementwise /~ win, and even more so for float[8]. (The first method in the float[8] set below is also less accurate, unless Digits is raised.)

restart:
with(LinearAlgebra): dtype:=anything:

a := RandomVector(500, generator = 1 .. 100, outputoptions=[datatype=dtype]):
b := RandomVector(500, generator = 1 .. 100, outputoptions=[datatype=dtype]):

st := time[real](): for j to 10^4 do a*~ (x -> 1/x)~(b) end do: time[real]()-st;
                             4.050

st := time[real](): for j to 10^4 do a*~ (1/~(b)) end do: time[real]()-st;
                             1.781

st := time[real](): for j to 10^4 do zip(`/`,a,b) end do: time[real]()-st;
                             1.567

st := time[real](): for j to 10^4 do a/~b end do: time[real]()-st;
                             1.358

div_vec := proc (a::Vector, b::Vector) 
local i, c;
 c := Vector(Dimension(a)): 
for i to Dimension(a) do c[i] := a[i]/b[i] end do: 
return c :
end proc:

st := time[real](): for j to 10^4 do div_vec(a, b) end do: time[real]()-st;
                             7.659

res1:=a*~ (x -> 1/x)~(b):
res2:=a*~ (1/~(b)):
res3:=zip(`/`,a,b):
res4:=div_vec(a, b):
res5:=a/~b:

Norm(res4-res1),Norm(res4-res2),Norm(res4-res3),Norm(res4-res5);
                           0, 0, 0, 0

restart:
with(LinearAlgebra): dtype:=float[8]:
a := RandomVector(500, generator = 1 .. 100, outputoptions=[datatype=dtype]):
b := RandomVector(500, generator = 1 .. 100, outputoptions=[datatype=dtype]):

st := time[real](): for j to 10^4 do a*~ (x -> 1/x)~(b) end do: time[real]()-st;
                             10.667

st := time[real](): for j to 10^4 do a*~ (1/~(b)) end do: time[real]()-st;
                             0.303

st := time[real](): for j to 10^4 do zip(`/`,a,b) end do: time[real]()-st;
                             0.347

st := time[real](): for j to 10^4 do a/~b end do: time[real]()-st;
                             0.176

div_vec := proc (a::Vector, b::Vector) 
local i, c;
 c := Vector(Dimension(a)): 
for i to Dimension(a) do c[i] := a[i]/b[i] end do: 
return c :
end proc:

st := time[real](): for j to 10^4 do div_vec(a, b) end do: time[real]()-st;
                             7.501

res1:=a*~ (x -> 1/x)~(b):
res2:=a*~ (1/~(b)):
res3:=zip(`/`,a,b):
res4:=div_vec(a, b):
res5:=a/~b:

Norm(res4-res1),Norm(res4-res2),Norm(res4-res3),Norm(res4-res5);
                        -9                        -15        
  6.66666721826913999 10  , 3.55271367880050093 10   , 0., 0.

You would re-use a procedure returned by dsolve/numeric by putting it into a Library archive, not a text file.

restart:

sol:=dsolve({diff(f(x),x)=sin(x),f(0)=1.0},
            numeric,output=listprocedure);\

            sol := [x = proc(x)  ...  end proc;, f(x) = proc(x)  ...  end proc;]

Y:=eval(f(x),sol):

Y(4.6);

                                    2.11215237701071

LibraryTools:-Create("c:/temp/abond.mla",'WRITE'):
libname:="c:/tmp/abond.mla",libname:

savelib(Y);

And now, in a new session

restart:

libname:="c:/tmp/abond.mla",libname:

Y(4.6);

                                    2.11215237701071

If you use the `known=F` option when computing `sol`, then you might have to savelib `F` as well (I haven't checked this).

Sorry that I haven't had any time to look at your followup question in an older message thread, about trying to speed up your earlier techniques or iterated solving.

You can use either the plots[fieldplot] or the DEtools[DEplot] commands.

This is a frequently asked questioning line.

There are several ways, none perfect. Choose which suits your more complicated examples.

See this reference, and follow its own links, for more ideas and methods.

For this simple example you posted, you could do something like this (with names as I have here, or with strings).

ans:=`2+2`;

                              2+2

parse(ans);

                               4

Below, I am deliberately forming the indexed `F` that you describe. But it's not necessary to do all that, for your example. The easiest alternative is to use the `plot3d` command directly. (There are various other ways, between these two.)

for i from 1 to 100 do
  x[i]:=(i-1)*0.1:
  y[i]:=(i-1)*0.2:
  F[i]:=x[i]^2 +x[i]*y[i]+y[i]^2;
end do:

plots:-surfdata( Array(1..100,1..100,(i,j)->F[i],datatype=float[8]),
                 0..9.9, 0..19.8, axes=box );

plot3d(x^2 + x*(2*x) + (2*x)^2, x=0..9.9, y=0..19.8, axes=box); 

for i from 1 to 100 do
 x[i]:=(i-1)*0.1:
end do:
for j from 1 to 50 do
 y[i]:=(i-1)*0.2:
end do:
for i from 1 to 100 do
 for j from 1 to 50 do
  F[i,j]:=x[i]^2 +x[i]*y[j]+y[j]^2;
 end do:
end do:

plots:-surfdata( Array(1..100,1..50,(i,j)->F[i,j],datatype=float[8]),
                 0..9.9, 0..9.8, axes=box );

plot3d(x^2 + x*y + y^2, x=0..9.9, y=0..9.8, axes=box);

> map(fnormal, E3);

                 [0.1745329247  0.6981317015  0.1745329256]

For your particular example, another way would be to not `evalf` the quantity assigned to `L`, but to leave as exact symbolic until the later.

> L := 2*Pi/(6*2):
 
> E3 := simplify(map(int, N0, eta = eta0 .. eta2));

                               [1      2     1    ]
                         E3 := [-- Pi  - Pi  -- Pi]
                               [18     9     18   ]
> evalf(E3);

                 [0.1745329252  0.6981317008  0.1745329252]