Preben's observation that those numeric integrals don't depend upon `t`, and can be computed just once up front (instead of repeatedly, for each step) is very good.

As a minor comment, that would also work with the original mode of using dsolve/numeric.

restart:
st:=time():
# Maple15\xsys11demo.mw, standard worksheet, BoltzQM, Shift & enter, ADD not SUM, 1-D math for procedure
with(StringTools):
FormatTime(%c);
with(DEtools):with(plots):
printlevel := 1:
mm := 3: nn := 3: a := 5.: e := 1.:  R := .1: h := 1.:

for n to nn do
for nd to nn do
for md to mm-1 do
integ[n,nd,md]:=evalf(Int( sin(Pi*n*x/a)*sin(Pi*nd*x/a)*sin(Pi*n*y/a)
                           *sin(Pi*nd*y/a)/sqrt((x-y)^2+4*R^2*sin(Pi*md/mm)^2),
                           y = 0. .. a, x = 0. .. a))
end do
end do
end do;

dproc3 := proc (nn, t, Y, YP) local n, nd, md;
  for n to nn do
    YP[n] := -((I*e^2*4)/(h*a^2))*Y[n]
             *add(add(Y[nd]*conjugate(Y[nd])*integ[n,nd,md],
                      md = 1 .. mm-1), nd = 1 .. nn);
  end do;
end proc:

ic3 := array([seq(1/sqrt(nn), n = 1 .. nn)]):
dvars3 := [seq(c[n](t), n = 1 .. nn)]:

p := dsolve(type=numeric, number = nn, procedure = dproc3, start = 0.,
            initial = ic3, procvars = dvars3, range=0..0.1, method=rkf45):

FormatTime(%c);
odeplot(p,[seq([t,abs(c[n](t))],n=1 .. nn)]);
time()-st;

@herclau

1) Your code does not produce the optimal value of Qnum in the plotted range, since the plot structure contains only value for only a discrete grid of points (`Do`,N`). Using Optimization:-Maximize gets a better value.

2) Your pasted code uses % in the last line, but you likely meant %%% (or else the posted code produces an error.) You probably intended to apply your procedure `maxindex` to the plot structure. But more seriously, it returns an index for the entry [14,5] which is incorrect and does not match the found maximal value of the plot structure. (That's probably because your code uses the only submatrix 2..,... the point of which isn't clear.)

p1 := plot3d(Qnum(Do,N),Do=0.05..0.5,N=1..800,axes=boxed,
             labels = ['Do','N', 'Q'], view = [0.05 .. 0.5,1..800, 0..25000]):

M := plottools:-getdata(p1)[3]:

max(%[2..,..]);

                      24532.2980300000017

maxindex := proc(A::Matrix) local i;
     op(attributes(max(seq(setattribute(Float(abs(rhs(i)), 0),
     [lhs(i)]), i = op(2, A))))) end proc:

maxindex(%%%[2..,..]); # gives the wrong index

                             14, 5

rtable_scanblock( M, [rtable_dims(M)],
    (val,ind,res) -> `if`(val > res[2],[ind,val],res),
    [[1,1],-infinity]);

                    [[15, 5], 24532.29803]

MM[15,5], MM[14,5];

                   24532.29803, 24216.38369

The Optimization:-Maximize result of 24560.24... was greater in value.

@herclau

1) Your code does not produce the optimal value of Qnum in the plotted range, since the plot structure contains only value for only a discrete grid of points (`Do`,N`). Using Optimization:-Maximize gets a better value.

2) Your pasted code uses % in the last line, but you likely meant %%% (or else the posted code produces an error.) You probably intended to apply your procedure `maxindex` to the plot structure. But more seriously, it returns an index for the entry [14,5] which is incorrect and does not match the found maximal value of the plot structure. (That's probably because your code uses the only submatrix 2..,... the point of which isn't clear.)

p1 := plot3d(Qnum(Do,N),Do=0.05..0.5,N=1..800,axes=boxed,
             labels = ['Do','N', 'Q'], view = [0.05 .. 0.5,1..800, 0..25000]):

M := plottools:-getdata(p1)[3]:

max(%[2..,..]);

                      24532.2980300000017

maxindex := proc(A::Matrix) local i;
     op(attributes(max(seq(setattribute(Float(abs(rhs(i)), 0),
     [lhs(i)]), i = op(2, A))))) end proc:

maxindex(%%%[2..,..]); # gives the wrong index

                             14, 5

rtable_scanblock( M, [rtable_dims(M)],
    (val,ind,res) -> `if`(val > res[2],[ind,val],res),
    [[1,1],-infinity]);

                    [[15, 5], 24532.29803]

MM[15,5], MM[14,5];

                   24532.29803, 24216.38369

The Optimization:-Maximize result of 24560.24... was greater in value.

@Alejandro Jakubi Yes, that is what I meant, about those last piecewise results (which I got on Windows 7, 64bit Maple 15.01, FWIW). Ie. Pi is not variable, and the nonzero piecewise part is disallowed by the assumptions on n.

I forgot to mention (and you probably already noticed)

> int(abs(cos(n*x)),x=0..Pi) assuming n::posint; 

                                       2

> int(abs(cos(n*x)),x=0..Pi) assuming n::negint;

                                       2

One thing that is curious about this example is that the only two methods that succeed are ftoc and ftocms, but the correct result can be obtained for them (using say n=13) for a nonnegative parameter b which is then subsequently evaluated at b=Pi.

Using Maple 15.01,

> restart:
> kernelopts(printbytes=false):

> ans:=int(abs(cos(13*x)),x=0..b) assuming b>0:

> hastype(ans,piecewise);                      

                                     true


> eval(ans,b=Pi);                              

                                       2


> int(abs(cos(13*x)),x=0..Pi);

                                      24
                                      --
                                      13

And in more detail,

> restart:

> int(abs(cos(13*x)),x=0..b,method=_RETURNVERBOSE) assuming b>0; # takes a minute

        4096                 12   11264                 10
[ftoc = ---- %1 sin(b) cos(b)   - ----- %1 sin(b) cos(b)
         13                        13

       11520                 8   5376                 6
     + ----- %1 sin(b) cos(b)  - ---- %1 sin(b) cos(b)
        13                        13

       1120                 4   84                 2
     + ---- %1 sin(b) cos(b)  - -- %1 sin(b) cos(b)  + 1/13 %1 sin(b)
        13                      13

       /{              -26 b + Pi                             \
       |{ 2/13 floor(- ----------) + 2/13        0 < 26 b - Pi|
     + |{                 2 Pi                                |, FAILS = (
       |{                                                     |
       \{                0                         otherwise  /

    distribution, piecewise, series, o, polynomial, ln, lookup, cook, ratpoly,

    elliptic, elliptictrig, meijergspecial, improper, asymptotic, meijerg,

    contour), ftocms = 1/13 sin(13 b) signum(cos(13 b))

       /{              -26 b + Pi                             \
       |{ 2/13 floor(- ----------) + 2/13        0 < 26 b - Pi|
     + |{                 2 Pi                                |]
       |{                                                     |
       \{                0                         otherwise  /

                        13               11               9              7
%1 := signum(4096 cos(b)   - 13312 cos(b)   + 16640 cos(b)  - 9984 cos(b)

                  5             3
     + 2912 cos(b)  - 364 cos(b)  + 13 cos(b))


> eval(%,b=Pi);

[ftoc = 2, FAILS = (distribution, piecewise, series, o, polynomial, ln, lookup,

    cook, ratpoly, elliptic, elliptictrig, meijergspecial, improper,

    asymptotic, meijerg, contour), ftocms = 2]


> restart:

> int(abs(cos(13*x)),x=0..Pi,method=_RETURNVERBOSE);

        24
[ftoc = --, FAILS = (distribution, piecewise, series, o, polynomial, ln,
        13

    lookup, cook, ratpoly, elliptic, elliptictrig, meijergspecial, improper,

                                            24
    asymptotic, meijerg, contour), ftocms = --]
                                            13

So, that's a bit of a discrepency, which suggests that the FTOC is not applied in that (naive, simple) manner.

I wonder whether there is anything special about using the constant name Pi. I tried raising Digits very high (10000) since `is` can rely upon `evalf`. But the wrong result reoccurred. These next results are a little curious, though, since `Pi` isn't a variable parameter:

> restart:

> int(abs(cos(n*x)),x=0..Pi,AllSolutions) assuming n::posint;  

                             /{                 Pi  \
                             |{ 1/n        Pi = --- |
                         2 + |{                 2 n |
                             |{                     |
                             \{  0         otherwise/

> restart:

> int(abs(cos(n*x)),x=0..Pi,AllSolutions) assuming n::negint;  

                            /{                     Pi \
                            |{ - 1/n        Pi = - ---|
                        2 + |{                     2 n|
                            |{                        |
                            \{   0          otherwise /

> restart:

> int(abs(cos(n*x)),x=0..Pi,AllSolutions) assuming n::{negint}; 

              /{                  Pi \   /{                     Pi \
              |{ 1/n        0 = - ---|   |{ - 1/n        Pi = - ---|
          2 - |{                  2 n| + |{                     2 n|
              |{                     |   |{                        |
              \{  0         otherwise/   \{   0          otherwise /

I only include that last one above because for n=13 the result of 24/13 is wrong by 2/n, and nothing else I've seen comes close to suggesting how that result is obtained. Of course, someone could also debug it and find out exactly what it does.

For the computation involving the norms, it would be better to first assign the adjusted Vectors trf(x) and trf(y) to variables, so that the computation is not repeated.

Using Marikyan's elementwise construction of the adjusted Vectors

with(Statistics): with(LinearAlgebra):
randomize():

y := Sample(Normal(0, 1), 10):
yadj:=y-~Mean(y):

x := Sample(Normal(0, 1), 10):
xadj:=x-~Mean(x):

cos(VectorAngle(xadj, yadj));

                                0.6246925138

Correlation(x, y);

                                0.6246925138

Norm(yadj,2)*cos(VectorAngle(xadj,yadj))/Norm(xadj, 2);

                                0.4077383202

LinearFit([1, t], x, y, t);

                  -0.155697529986061 + 0.407738320274666 t

This is also an interesting example for illustrating the relative cost of loops versus in-place operations versus elementwise operations on float[8] rtables.

with(Statistics): with(LinearAlgebra):
randomize():

Digits:=15:

N:=10^6:

y := Sample(Normal(0, 1), N):
yadj:=y-~Mean(y):
x := Sample(Normal(0, 1), N):
xadj:=x-~Mean(x):

cos(VectorAngle(xadj, yadj));
                             0.00120391456372874

Correlation(x, y);
                             0.00120391456372511

Norm(yadj,2)*cos(VectorAngle(xadj,yadj))/Norm(xadj, 2);
                             0.00120342424655244

LinearFit([1, t], x, y, t);
                -0.000845803685442095 + 0.00120342424654874 t

trf := proc(a::Vector)
  # transformation of a vector to a new vector
  # centred about it's mean, with length equal
  # to the original vector's standard deviation
  uses LinearAlgebra, Statistics:
  local i,b,m,n:
  m:=Mean(a): n:=Dimension(a): 
  b:= Vector[row](n):
  for i from 1 to n do
    b[i]:= evalf((a[i]-m)/(sqrt(n-1))):
  end do:
  return b:
end proc:

res:=CodeTools:-Usage(trf(x)):
memory used=2.44GiB, alloc change=0 bytes, cpu time=45.43s, real time=45.47s

Norm(res-xadj/evalhf(sqrt(N-1)));
                                                -17
                          2.86229373536173171 10   

trf := proc(a::Vector)
  uses LinearAlgebra, Statistics:
  local i,b,m,n,sn:
  m:=Mean(a): n:=Dimension(a):
  sn := evalhf(sqrt(n-1));
  b:= Vector[row](n,i->(a[i]-m)/sn,datatype=float[8]):
  return b;
  end proc:
res:=CodeTools:-Usage(trf(x)):

memory used=76.31MiB, alloc change=0 bytes, cpu time=1.64s, real time=1.64s

Norm(res-xadj/evalhf(sqrt(N-1)));
                                                -17
                          2.86229373536173171 10   

trf := proc(a::Vector)
  uses LinearAlgebra, Statistics:
  local i,b,m,n:
  m:=Mean(a): n:=Dimension(a):
  b:= copy(a):
  map[evalhf,inplace](`-`,b,m);
  VectorScalarMultiply(b,evalhf(1/sqrt(n-1)),inplace);
  return b;
  end proc:

res:=CodeTools:-Usage(trf(x)):
memory used=7.64MiB, alloc change=0 bytes, cpu time=31.00ms, real time=31.00ms

Norm(res-xadj/evalhf(sqrt(N-1)));
                                                -17
                          2.86229373536173171 10   

CodeTools:-Usage(x-~Mean(x)):
memory used=7.64MiB, alloc change=0 bytes, cpu time=15.00ms, real time=15.00ms

CodeTools:-Usage(xadj/evalhf(sqrt(N-1))):
memory used=7.63MiB, alloc change=0 bytes, cpu time=0ns, real time=0ns

Adding a view=0..1000 option to the plots:-display(...) calls at the end of the "Action When Value Changes" code of all four sliders might make for a a little easier visual insight.

Also, while it does not make much difference for this small and naturally quickly solved example, for larger and more computationally intensive problems there can be a benefit from using odeplot and fieldplot instead of DEplot.

Adding a view=0..1000 option to the plots:-display(...) calls at the end of the "Action When Value Changes" code of all four sliders might make for a a little easier visual insight.

Also, while it does not make much difference for this small and naturally quickly solved example, for larger and more computationally intensive problems there can be a benefit from using odeplot and fieldplot instead of DEplot.

Evalf at `Digits` worth of digits of working precision does not promise any specific number of digits of accuracy for floating-point evaluation of a compound expression.

On the other hand, is n always an integer?. Exact simplification might be possible.

Perhaps you could paste the expression in a 1D text (after applying the lprint command, or upload it in a worksheet). That saves people from each having to type in your formula from the embedded image.

@amrramadaneg As I mentioned, there are various ways to get the same effect. It's still not clear in what format your data is given. (Also, you began with 3D data, but now you talk of `v` and calls to the 2D `plot` command.)

Another couple of ways to do the second example:

restart:
x:=Vector(100,i->(i-1)*0.1):
y:=Vector(50,i->(i-1)*0.2):

F:=Array(1..100,1..50,(i,j)->x[i]^2 +x[i]*y[j]+y[j]^2,datatype=float[8]):

plots:-surfdata(F,0..9.9,0..9.8,axes=box);

A:=Array(1..100,1..50,1..3,(i,j,k)->`if`(k=1,x[i],`if`(k=2,y[j],F[i,j])),datatype=float[8]):
plots:-display(PLOT3D(MESH(0..9.9,0..9.8,0..600,A)),axes=box);

@amrramadaneg As I mentioned, there are various ways to get the same effect. It's still not clear in what format your data is given. (Also, you began with 3D data, but now you talk of `v` and calls to the 2D `plot` command.)

Another couple of ways to do the second example:

restart:
x:=Vector(100,i->(i-1)*0.1):
y:=Vector(50,i->(i-1)*0.2):

F:=Array(1..100,1..50,(i,j)->x[i]^2 +x[i]*y[j]+y[j]^2,datatype=float[8]):

plots:-surfdata(F,0..9.9,0..9.8,axes=box);

A:=Array(1..100,1..50,1..3,(i,j,k)->`if`(k=1,x[i],`if`(k=2,y[j],F[i,j])),datatype=float[8]):
plots:-display(PLOT3D(MESH(0..9.9,0..9.8,0..600,A)),axes=box);

Of course, you could also get a similar discretization from plot3d and its `grid` option.

plot3d(x^2 + x*y + y^2, x=0..9.9, y=0..9.8, grid=[100,50], axes=box);

Of course, you could also get a similar discretization from plot3d and its `grid` option.

plot3d(x^2 + x*y + y^2, x=0..9.9, y=0..9.8, grid=[100,50], axes=box);

@kutscher 

fsolve(ddEQ2,0.0..1.0);

                    0.089572502098965952134

plot(ddEQ2,0.0..0.3,view=-1000..1000);

plot(ddEQ2,0.3..1.0,view=-1..1);