Question: Thanks, and how can I make and use functions f( c[n](t) ) of the basic output c[n](t) of dsolve/numeric, e.g. norm = add(abs(c[n](t))^2,n=1..nn) ?

# Maple15\xsys14demo.mw, Maplesoft, standard worksheet, BoltzQM, Shift & enter, ADD not SUM, 1-D math for procedure

restart;
with(StringTools):
FormatTime(%c);
with(DEtools):
with(plots):
printlevel := 1:
mm := 3:
nn := 3:
a := 5.:
e := 1.:
R := .1:
h := 1.:
t2:=5:
Δt:=t2/10:

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])
          *evalf[11](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])), md = 1 .. mm-1), nd = 1 .. nn);
   end do;
end proc;
ic3 := array([seq(((1-I)/sqrt(2))/sqrt(nn), n = 1 .. nn)]):
dvars3 := [seq(c[n](t), n = 1 .. nn)]:
tt := time():
p := dsolve(type=numeric, number = nn, procedure = dproc3, start = 0., initial = ic3, procvars = dvars3,
     range=0..t2+Δt, method=rkf45);

p(3);
p(3)[2];
p[2](3);

evalf[11](add(abs(c[n](3))^2,n=1..nn)); (* norm, not evaluated error *)


#abs(evalf[11](p(3)[2]));
#Int(2.*ln(abs(p[3](t))), t = 5. .. 5.5000000000);

`Solver time = `,time()-tt;
FormatTime(%c);
tt := time():
display([seq(odeplot(p,[t,ln(abs(c[n](t))^2)]),n=1 .. nn)]);


display([seq(odeplot(p,[t,evalf[11](Int(ln(abs(c[n](td))^2),td= t..t+Δt))]),n=1 .. nn)]); (* running average, error message*)


evalf[11](Int(ln(abs(c[3](t))^2),t= t2..t2+Δt));
`Plotting time = `,time()-tt;

Error, (in plots/odeplot) curve is not fully specified in terms of the ODE solution, found additional unknowns {td, c[1](td)}