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

If both the following redefinitions are made to `plots/animate` then animate works:

`plots/animate`:=subs(subs=((x,y)->eval(y,x)),eval(`plots/animate`)):
`plots/animate`:=subs(evalf=(()->args),eval(`plots/animate`)):

animate(PointPlot, [L[t]], t=1..5,frames=5,digits=13);
animate(PointPlot, [L[round(t)]], t=1..5,paraminfo=false);

While the first change is safe and works fine (in my opinion - I have used it a lot since 2007), the last one surely is not, but it points to the problem.

In this case the last change makes the title weird looking. I have set digits=13 just to call attention to the title problem That issue is, however, easily fixed.

Preben Alsholm

restart;
A:=Matrix(3,5,symbol=a);
<seq(LinearAlgebra:-Row(A,i),i=combinat:-randperm(3))>;

Preben Alsholm

How about putting all of your procedures into a module:

M:=module() export a,b,c,..,d; local A,B,C, ..., D; global g1, g2, ..., gn;
a:=proc(x,y) ....end proc;
b:=proc(x,y,z) ... end proc;
....
end module;

Preben Alsholm
 

If you use piecewise, life is easier.

restart;
eq:=diff(q(t), t) =piecewise(t<8760,3-3*q(t)/10^4, -3*q(t)/10^4);
#First version:
p1:=dsolve({eq, q(0)=0}, type=numeric, range=0..16000):
plots[odeplot](p1,refine=1);
#Second version:
DEtools[DEplot](eq,q(t),t=0..16000,[q(0)=0],linecolor=blue,color=gray);

Preben Alsholm
 

Doug Meade's version may be written somewhat shorter this way:

XY := proc(T,H)
  global integ;
  if not ( T::numeric and H::numeric ) then return ['procname'(args)[1],'procname'(args)[2]] end if;
  integ(parameters=[H]);
  eval( [x(t),y(t)], integ(T) );
end proc:
XY(2,1);
                 [2.08917117586676548, 1.35878133507557242]
XY(u,8);
                         [XY(u, 8)[1], XY(u, 8)[2]]

plot3d( XY(t,H), t=0..3, H=-1..5, color=[red,blue], transparency=0.5, axes=boxed );

It should be noticed that when XY returns 'unevaluated', it returns a list, just as in the case of numeric input. (This means that strictly speaking, it doesn't return unevaluated).

Otherwise the plot3d won't work.

Simple tests (see below) seems to indicate that the XY- version takes about half the time compared to the separate X, Y - version.
Not really surprising.

Preben Alsholm

t0:=time():
for k to 100 do
   for j to 100 do
    XY(j*.1,k*.1)
   end do
end do:
time()-t0;
t0:=time():
for k to 100 do
   for j to 100 do
    [X(j*.1,k*.1),Y(j*.1,k*.1)]
   end do
end do:
time()-t0;

 

If you are using 2D-math, then the following assignment

f(t):= t^2

is considered ambiguous by Standard Maple.

If you intend to define a function, then the proper definition is

f:= t -> t^2;

If you use Maple Notation for your input (which I always do, and ask my students to do), then the assignment

f(t):= t^2;

is writing to the remember table for the function f. This means that f now knows what to do if it literally receives the letter t. It won't know what to do about anything else, unless of course you have already defined f as a function (= procedure), as in

restart;

f:= t -> t^2;

f(t):=411;

Now try

f(t);

f(6),

f(x);

Besides asking my students to use Maple input, I ask them to use Worksheet interface.

These setting can be changed by going to  Tools/Options/Display  and Tools/Options/Interface.

Preben Alsholm

If you use the long form of polarplot, plots:-polarplot (or plots[polarplot]) the plot is fine:

Explore(plots:-polarplot(1+a*cos(theta), theta = 0 .. 2*Pi)) ;

Preben Alsholm

I'm not sure what you want to do, but is this anywhere close?

for k to 10 do fnormal(10^(-k),7) end do;
 

Preben Alsholm

Another (obvious) solution is to place with(plots): outside the loops.

Preben Alsholm

I have changed the syntax for the initial and boundary conditions and used Pi instead of pi (in this case not important as you will notice, since I replace it by evalf(Pi)).

Also all u's have been made u[0]'s.

eq_u_0 := diff(u[0](x,t),`$`(x,4))+a*u[0](x,t)+diff(u[0](x,t),`$`(t,2))=0;
bc_u0 :=u[0](0,t) = 0, u[0](Pi,t) = 0, D[1,1](u[0])(0,t) = 0, D[1,1](u[0])(Pi,t) = 0, u[0](x,0) =a*sin(x), D[2](u[0])(x,0) = 0;
#Exact solution, no conditions:
pdsolve(eq_u_0);
#Numerical solution (for a = 1):
M:=pdsolve(op(eval([eq_u_0,{bc_u0}],{a=1,Pi=evalf(Pi)})),numeric);
M:-plot(x=0..Pi,t=1);
M:-animate(x=0..Pi,t=0..10,frames=100);
 

Preben Alsholm

I added two options (in Maple 13.02):

ff:=dsolve(syst union CI,{x(t),y(t)},type=numeric,output=listprocedure,maxfun=0,range=0..0.1);

fx := subs(ff,x(t)): fy := subs(ff,y(t)):
fx(0), fy(0);
fx(.01);
plot(fx,0..0.1);
plot(fy,0..0.1);
 

The plots are rather wild.

Preben Alsholm