At your level of Maple, I think it's a bad idea to try writing functions such as those you quoted, that are very much specialized for a particular data structure.  It's better to look at the problem at a more mathematical, rather than structural, level.  As I understand it, you have an expression in sines and cosines with one variable t, perhaps something such as this

> S := 3*sin(2*t) + sin(3*t + 2) + cos(t)^2;

The first step is to take the Fourier transform.  This can be done with fourier in the inttrans package.

> with(inttrans):

> F := fourier(S, t, w);

F := 1/2*Pi*(Dirac(w-2)-6*I*Dirac(w-2)+2*Dirac(w)+Dirac(w+2) +6*I*Dirac(w+2)+2*I*exp(-2*I)*Dirac(w+3)-2*I*exp(2*I)*Dirac(w-3))

Next we collect terms involving the same Dirac function.

> F := collect(F, Dirac);

F := (1/2-3*I)*Pi*Dirac(w-2)+(1/2+3*I)*Pi*Dirac(w+2)+ exp(-2*I)*Pi*Dirac(w+3)*I-I*exp(2*I)*Pi*Dirac(w-3)+Pi*Dirac(w)

We make this into a list.  This step wouldn't work if there was just one term, but I won't worry about that.

> L := convert(%, list);

L := [(1/2-3*I)*Pi*Dirac(w-2), (1/2+3*I)*Pi*Dirac(w+2), exp(-2*I)*Pi*Dirac(w+3)*I, -I*exp(2*I)*Pi*Dirac(w-3), Pi*Dirac(w)]

Each item in the list is of the form c * Dirac(w - p), and I want to make that into a list [p,w] representing a point to plot.  I'll define a procedure that does that to an item.  Extracting the c and p will use Maple's pattern-matching capabilities.

> GetPoint := proc(Term)  
      if match(Term = c*Dirac(w-p), w, 's') 
      then eval([p, abs(c)], s)
      end if
  end proc;

We use map to apply this to every item in the list L.

> Points := map(GetPoint, L);

Points := [[2, 1/2*37^(1/2)*Pi], [-2, 1/2*37^(1/2)*Pi], [-3, Pi], [3, Pi], [0, Pi]]

And now the result can be plotted.

> plot(Points, style = point);

Although I'm not familiar with MADONNA, this sounds to me like a DAE (differential-algebraic equations) problem. I think dsolve(...,numeric) should be able to handle it.  See the help page ?dsolve,numeric,DAE.

diff(y(t), t, t) = piecewise(y(t)+diff(y(t), t) > 0, 1, y(t)+diff(y(t), t) < 0, -1)

You might also try symbolic solution of your DE.

dsolve can do nothing with it, because of the piecewise.   You need to look at the two pieces separately.

Suppose e.g. you start at y(0) = -2.1 y'(0) = 2.  Thus y(0) + y'(0) = < 0, so this is in the second part.

sol1 :=  dsolve([diff(y(t),t,t)=-1,y(0)=-2.1,D(y)(0)=2]);

sol1:=y(t)=-1/2*t^2+2*t-21/10

Now what does that make of y + y' ?

q1 := eval(y(t)+diff(y(t),t), sol1);

q1:=(1/2)*t^2+t-1/10

And when does that hit 0?

solve(q1 = 0);

1-2/5*sqrt(5),1+2/5*sqrt(5)

The first of these is encountered first as t increases from 0.

What happens outside this interval?  Well, strictly speaking, your DE is not defined when y + y' = 0.  But we can try to match sol1 with a solution of the other piece with the same values of y and y' at the boundary point.

t0:= %[1];

sol2:= dsolve([diff(y(t),t,t)=1,y(t0)=eval(rhs(sol1),t=t0),  
    D(y)(t0)=eval(diff(rhs(sol1),t),t=t0)]);

sol2 := y(t)= (1/2)*t^2+(4/5)*sqrt(5)*t-3/10-(4/5)*sqrt(5)

What does that make of y + y'?

q2 := eval(y(t)+diff(y(t),t), sol2);

q2:=(1/2)*t^2+(4/5)*sqrt(5)*t-3/10+t

This turns out to be positive for all t > t0, so sol2 gives the solution there.

The short-cut syntax for Vectors is with angle brackets < >.  For example:

> P := < 1, 2>;

You can then add Vectors using +, or multiply by scalars using *: for example,

> 3*P - <1,-1>;

The problem with the plot in Classic is that your y values are too large for the plot renderer.  You might try a re-scaling, e.g. use k!/N! instead of k!.

int(x^n*exp(-x),x=0..infinity) is n! when n is a nonnegative integer, but the integral exists for all n >= 0 (in fact for complex numbers with Re(n) > 0).  Thus the value using the Gamma function is more general.  You can use convert(%, factorial).

For the 3d plot: assume(n::posint) does not affect plot3d, so if you want only positive integer values of n you'll need to do something else.  Also note that if you define f as a function, you can either use

plot3d(f(x,n),  x = 0 .. 10, n = 1 .. 10)

or

plot3d(f, 0 .. 10, 1 .. 10)

It's not too hard to program the Newton-Raphson method yourself.  Something like this, perhaps.  This assumes it is possible to take derivatives of your procedures.
Unfortunately, JACOBIAN doesn't work properly for your example using fsolve.  You could perhaps write your own code to produce the Jacobian matrix.
 

J:= Matrix @ codegen[JACOBIAN]([f, g]);
Newt := X -> X - J(X[1],X[2])^(-1) . <f(X[1],X[2]), g(X[1],X[2])>;
X[0] := <5.0, 6.0>;
for count from 1 to 5 do
  X[count]:= Newt(X[count-1])
end do;

I suppose you mean "object" in the sense of object-oriented programming.  Yes, this can be done.  Start by looking at the help page ?record.

Assuming it's the sum of the squares of the distances from the data points P_j to the line that you wish to minimize, the solution can be described parametrically as the line X = A + t B where B is a nonzero vector.  If I'm not mistaken, we can take A to be the average of the P_j, and then B should be an eigenvector for the largest eigenvalue of the matrix M = sum_j (P_j - A) (P_j - A)^T (where the vectors are taken to be column vectors).  For example:

with(LinearAlgebra):
Data:= Matrix([[1.33, 2.57, 2.24], [3.68, 5.27, 6.99], 
[0.92, -0.98, 0.27], [-0.72, -0.74, 0.8], [0.29, -0.31, -0.77], 
[-0.93, -0.2, -0.4], [0.69, 0.44, 0.67]]);
n:= RowDimension(Data);
A:= Data^%T . < (1/n $ n) >;
Q := Data^%T - A . < (1 $ n)>^%T ;
M := Matrix(Q . Q^%T, shape=symmetric);
E,V:= Eigenvectors(M);
B:= Column(V, 3);

And to visualize it.

with(plots):
dataplot:= pointplot3d(Data, symbol=box,colour=black):

tmin,tmax:= (min,max)(seq(B^%T . Column(Q,j), j=1..n));
lineplot:= plottools[line](convert(A+tmin*B,list),  
   convert(A+tmax*B,list), colour=red):
display([dataplot,lineplot], scaling=constrained);

I'm a bit confused about what you want to do here.  Do you want the procedure to pretty-print 2D math, something like this?

> p1:= x -> '`` implies (F = x)';

Or do you want expressions in 1D text, like this?

> p2 := x -> printf("=> F = %a\n", x);

The fsolve that is causing the error message is the one in the procedure g: unfortunately your g(x,y) calls the procedure g with symbolic x and y, and that won't work.  fsolve could deal with a single procedure; unfortunately it can't yet deal with more than one procedure. 

fsolve({f,g},[-10..10,-10..10]);

  Error, (in fsolve) Case of systems of equations involving procedures is not implemented yet
 

Your insistence on using the first solution computed by fsolve (do you really mean that?) makes things a bit difficult: if you were willing to have g(x,y) = 0 if _any_ root of q^5 + x*q^2 + y*q had real part -1, you could just say

geqs:= evalc({Re,Im}( (-1+I*t)^5 + x*(-1+I*t)^2 + y*(-1+I*t)));

fsolve({f(x,y)=0,op(geqs)},{x,y,t},{x=-10..10,y=-10..10});

{x = -.9727007669, y = -1.972700767, t = 0.}

For next time, it would be better to post a worksheet rather than a .pdf file, so we could have the actual commands to work with without doing a lot of copy-and-paste.

I think there's a bug in the DAE solver here: if you execute the same command

f1 := dsolve(...);

twice, you get a solution the first time but an error the second time.

It seems to be related to the presence of the relerr option.  I'm filling out an SCR for this.

That does seem to be a limitation of the plot renderer in Maple 10 (Standard), but not in Maple 11.  Compare the following plots, both results of

plot(exp(x), x=-20..0, axis[2]=[mode=log]);

(Maple 10)

(Maple 11)

The actual plot structures are essentially identical.

Expanding on John May's code:

ftmp:=convert(f,pwlist);
flist:=[seq](ftmp[i], i=1..nops(ftmp), 2);
boundlist:= [-infinity,seq(ftmp[i],i=2..nops(ftmp),2),infinity];
critpoints:= [seq](fsolve(diff(flist[i],x),
                          x=boundlist[i]..boundlist[i+1]),
                               i=1..nops(flist));

Note that this gives you a list of critical points, some of which may be points of inflection rather than local minima or maxima; also, if the spline is not differentiable, some of the interval boundaries may be local minima or maxima.