The next problem is in your initial/boundary conditions.  You can't specify a condition at a single point, e.g. (L,0), it would have to be at (L,t) or (x,0).

You can't do it directly.  If your system involves f(anything(t)) as well as f(t) itself, you have functional differential equations rather than ordinary differential equations.  This makes a big difference, both theoretically (existence and uniqueness of solutions may be hard questions) and practically (dsolve won't touch it at all).  Perhaps if you posted your actual system somebody might have some idea of how to deal with it.

Hint for (a): The general solution of the DE is

where the sum is over all 2011'th roots of -1.  For this to be real,


For t -> infinity, the solution is dominated by the nonzero terms where Re(omega) is maximized.  These have frequency given by Im(omega)/(2*Pi).  For the largest possible frequency, take the roots closest to I and -I...

Yes, it's a bug.  As a work-around, you might try avoiding the use of the parameters option, or else use a different method.  If you're using gear because it's a stiff system, you might try  method=lsode.

Your code appears to be finding the maximum and minimum in the Vector y (which is not the same as local maxima and minima).  That can be done using rtable_scanblock as follows:

> rtable_scanblock( y, [],
     (val,ind,res) -> `if`(val > res[2],[ind,val],res), [[1],y[1]],
     (val,ind,res) -> `if`(val < res[2],[ind,val],res), [[1],y[1]]
 );

> plots[display](plottools[scale]~({p1,p2,p3},1/(2*Pi*sqrt(2)),1),
            labels=[x/(2*Pi*'sqrt(2)'),h]);

It can be done using subs.  Thus:

> f := x -> x*t;
   f3:= subs('t' = 3, eval(f));
   t:= 4;
   f3(x);

The quotes are only needed if the variable t has been assigned a value when f3 is being defined.

csgn is a complex version of the signum function.  It is 1 for a number in the right half-plane and -1 for a number in the left half-plane.   It will typically occur when using evalc on an expression containing square roots of complex expressions.  For example:

> evalc(sqrt(a+b*I));

Note that any expression for this will have to be discontinuous on the negative real axis (which is the branch cut of the principal branch of the square root), and this discontinuity is provided by the csgn.

In your case, the result of evalc contains

and

Assumptions such as beta > 0, Q > 0 should remove the need for csgn.

Perhaps something like this.

> n:= 2^10;
   M:= Array(0..n-1,0..n-1,(i,j) -> `if`((i-n/2)^2 + (j-n/2)^2 < 100,1,0),
     datatype=float[8]);
   F:= DiscreteTransforms:-FourierTransform(M);
   Fa:= abs~(F);
   Fs:= Array(1..n,1..n,datatype=float[8]);
   for i from 0 to 1 do for j from 0 to 1 do
      Fs[1+n/2*i .. n/2 + n/2*i, 1+n/2*j .. n/2 + n/2*j]:= 
      Fa[1+n/2*(1-i) .. n/2 + n/2*(1-i), 1+n/2*(1-j) .. n/2 + n/2*(1-j)]
   end do end do;
   Img:= ImageTools[Create](Fs,fit=true);
   ImageTools[View](Img);

If I understand the Matlab program, you want something like this:

> step:= 2*Pi/2000; 
   start_point1:=Pi/3+Pi/6;
   duration1:=Pi/3;
   start_point2:=4*Pi/3-Pi/6;
   duration2:=Pi/3;
   f:= t -> piecewise(t >= start_point1 and t <= start_point1 + duration1,2,
                  t >= start_point2 and t <= start_point2 + duration2,1,
                  0);
   plot(f(t),t=0..2*Pi);
   data:= Array(0..2000, i -> f(i*step));
   Y:= DiscreteTransforms:-FourierTransform(data);
   plot(<seq(<(i-1)*2*Pi/2000| abs(Y[i])>, i=1..2001)>);  

      

If I understand what you're trying to do:

> with(plots):
   n:= 2^10;
   implicitplot((X-n/2)^2+(Y-n/2)^2>=10^2,X=1..n,Y=1..n,grid=[100,100],
      gridrefine=2,filledregions=true, coloring=[black,white],title="Circular Aperture");

As Axel said, polynomials are probably the wrong class of functions.  I'll try constant + linear + exponential.

> fit2:= Fit(a + b*x + c*exp(-d*x),X,Y,x,parameterranges=[d=0..0.1]);

Of course you can't "deduce" it if you only have a finite number of terms, without knowing what the rest of the terms will be, but guessgf in the gfun package can sometimes make a good guess.  For example:

> with(gfun):
   S:=series(1-t+t^2-t^3,t,4);
   L:= seriestolist(S);
   guessgf(L,t,['ogf']);

You did provide chrem with two lists, namely [y,yp] and [1,3].  But the first one was a list of Vectors.  The type-checking does not check the type of the elements of the first argument; evidently a list of Vectors does not work for the first argument, though a list of lists would be ok.
Specifically, it seems the problem occurs because your first modulus is 1, and the code has a line

if m = 1 then v := 0 end if;

Later in the code, v will be subtracted from u which in your case will be a Vector. Addition or
subtraction of a Vector and 0 does not work, although with a List and 0 it would.

> q:= 0:
q + [1,2];

>  q + < 1,2 >;

Error, (in rtable/Sum) invalid arguments

(I'm assigning the value 0 to the variable q here, otherwise automatic simplification would 
change 0 +  to )

In this problem it suffices to consider f on the additive group generated by 1/20.  Unfortunately rsolve doesn't work for difference equations where the increment is a fraction, but we can get around that.

> eq:= f(x+y) = f(x) + f(y) + 80*x*y;
   geq:= eval(eq, f = (t -> g(20*t)));
   g1eq:= subs(x=t/20,y=1/20,g(1)=g1,geq);

Solve this equation for g:

> gt:=unapply(normal(rsolve({g1eq,g(1)=g1},g(t))),t);

Check that this satisfies the full equation

> normal(eval((lhs-rhs)(geq), g=gt));

Also check that g(1) = g1.

> gt(1);

Determine g1 from the condition f(1/4) = g(5) = 2.

> g1:= solve(gt(5) = 2, g1);

And now the answer:

> gt(16);