That was discussed here a few weeks ago: www.mapleprimes.com/forum/rangeanddomain

Maybe you don't need FunctionChart at all.  Something like this, perhaps?

> with(Statistics): with(plots):
    m, s := 6.8,1.47:
    f:= unapply(PDF(Normal(m,s),x),x):
    p1:= plot(f(x), x=0..m-s, filled=true,colour=orange):
    p2:= plot(f(x), x=m-s .. m+s, filled=true,colour=blue):
    p3:= plot(f(x), x=m+s..12, filled=true,colour=orange):
    p4:= pointplot({[m-s,f(m-s)],[m+s,f(m+s)]}, symbolsize=15):
    p5:= textplot([m-s-.3,f(m-s),mu-sigma],align=left),
        textplot([m+s+.3,f(m+s),mu+sigma],align=right):
  display([p1,p2,p3,p4,p5], xtickmarks=[m-s=mu-sigma,m=mu,m+s=mu+sigma]);

There seems to be a bug in Student[Calculus1][InflectionPoints] that prevents it from working with expressions involving floats, unless the option numeric is used.  Try

> Student[Calculus1]:-InflectionPoints(PDF(N,x), x=1 .. 12.5, numeric);

As for the graph: do you mean something like this?

>Student[Calculus1]:-FunctionChart(PDF(N,x), 1 .. 12.5, 
    xtickmarks=[%[1]=mu-sigma, 6.8=mu, %[2]=mu+sigma]);

1) Basically, you need to solve for as many variables as you have equations.  So: 

> S:= solve(q, {ds, k, m, n});

2) As you say, the result is a set.  The simplest way to get the result from a set is to use eval.  Thus:

> eval(ds, S);

3) Well, in this case you could try

> solve(q, i);

piecewise(j <= 2*m-4,[{i = 1/2*j+1}],[])

An animation in Maple has a fixed set of frames; if you run it in loop mode, it will repeat unchanged.  Perhaps what you want is something like this Maplet:

> with(Statistics) : with(Maplets[Elements]): randomize() :
  X:= RandomVariable(Normal(0,5)):
  makeanim:= proc(n)
   uses plots;
   local r,s,i,p,An;
   r:= Sample(X,n);
   s[1]:= 100;
   for i from 2 to n do s[i]:= s[i-1]+r[i] end do;
   for i to n do p[i]:= [seq(s[i],i=1..i)] end do;
   display([seq(LineChart(p[i],labels=["Time","Stock Price"],title=('time'=i)),i=1..n)],insequence=true);
  end proc: 
 M:= Maplet([Plotter[P](makeanim(50),play=true,continuous=false),
  [Button("Again",Evaluate('P'='makeanim(50)')),
  Button("Done",Shutdown())]]):

Maplets[Display](M);

To get the real part of each element of the Matrix M, you can try

> map(evalc @ Re, M);

I'm not sure what's producing your error message: perhaps you might post the actual code.  But it seems Eigenvalues can't handle a Matrix whose entries are polynomials involving complex numbers and floats.  For example:

> with(LinearAlgebra):
   M:= <<0.1, 0>|<0, I*x>>;
   Eigenvalues(M);

Error, (in LinearAlgebra:-LA_Main:-Eigenvalues) expecting either Matrices of rationals, rational functions, radical functions, algebraic numbers, or algebraic functions, or Matrices of complex(numeric) values

I sent in an SCR.

You can get around this (in a small case such as this, where numerical problems shouldn't be significant) by

> solve(CharacteristicPolynomial(M,r), r);

Something like this.  Note that I is used for the imaginary unit, so I'll change that to J.  It's also bad to have both a beta and a beta[T], so I'll just change the beta manually to 1.  You don't say what K and alpha are, so I'll take alpha as 1 and K as 0.08.  You also need an initial condition, which I'll choose arbitrarily.

> sys:= eval({diff(A(t),t) = phi*(1-A(t)/K)*F(t)-(gamma+mu[A])*A(t),
                        diff(J(t),t) = r*gamma*A(t) - M(t)*J(t)/(M(t) + M[T](t)) - 
                             beta[T]*M[T](t)*J(t)/(M(t)+M[T](t)) - mu[I]*J(t),
                        diff(F(t),t) = M(t)*J(t)/(M(t)+M[T](t)) - mu[F]*F(t),
                        diff(M(t),t) = (1-r)*gamma*A(t) - mu[M]*M(t),
                        diff(M[T](t),t) = alpha - mu[T]*M[T](t)},
               {phi = 5, mu[A]=0.05, mu[I]=0.05, mu[F] = 0.05, mu[U] = 0.05, mu[M]=0.1, 
                mu[T]=0.1, gamma = 0.075, r = 0.5, beta[T] = 0.3, alpha = 1, K = 0.08});

> ics:= {A(0)=1,F(0)=1,J(0)=2,M(0)=1,M[T](0)=0};
   S:= dsolve(sys union ics, numeric);

 > with(plots):
   odeplot(S,[[t,A(t)],[t,F(t)],[t,J(t)],[t,M(t)],[t,M[T](t)]], t=0..2,axes=box,
      colour=[red,green,blue,cyan,black]);

The problem is this:

sum(sum((-1)^(k-j)*binomial(k,j)*j^n, k=1..n ), j=1..k)

is using k in two different ways.  The k in k=1..n is a "dummy variable".  The k in j = 1 .. k is apparently a global variable.

It's actually surprising, and could be considered a bug, that when n is a positive integer Maple 12 gets a numerical value for this sum that
doesn't depend on k.  The reason, it seems, is that when j is a positive integer

sum((-1)^(k-j)*binomial(k,j)*j^n, k=1..n ) = 0 for j > n

Maple 12 evaluates the sum for j from 1 to k as the sum for j from 1 to n, which is also the sum for j from 1 to m for any m >= n.  I don't know if this is what catalan intended: I rather suspect that catalan made a mistake, and either both sums were supposed to be from 1 to n or the sum for j from 1 to k was supposed to be on the inside and the sum for k from 1 to n on the outside. 

Also make sure that you put that exp(x)*(x+1) in the denominator: you want
(-exp(-1)-x*exp(x))/(x*exp(x)*(x+1)), not (-exp(-1)-x*exp(x))/x*exp(x)*(x+1)
which would only have x in the denominator.
 

Maple doesn't support pseudographs, but it does support weighted graphs.  You could use edge weights to correspond to multiple edges, and vertex attributes to correspond to loops.  So something like this:
 

> G:= Graph([a,b,c,d,e],
{[{a,b},3],{a,e},{b,d},{b,e},{b,c},[{c,d},3],
{d,e}});

> SetVertexAttribute(G, a, [loops=1]):
  SetVertexAttribute(G, c, [loops=1]):
  for i in Vertices(G) do 
     if GetVertexAttribute(G,i,loops)=FAIL then
        SetVertexAttribute(G,i,[loops=0]) 
  end if end do:

> printf("This is a pseudograph with %d vertices and %d edges",
       nops(Vertices(G)),
       add(GetEdgeWeight(G,p), p=Edges(G)) 
        + add(GetVertexAttribute(G,v,loops), v = Vertices(G)));

> for v in Vertices(G) do
     printf("The degree of %a is %d\n", v,
      add(GetEdgeWeight(G,p), p=IncidentEdges(G,v)) +    
            2*GetVertexAttribute(G,v,loops))
  end do:

Actually the solve command works, but its result is three enormously complicated expressions.

>S:= [solve(%,rho)]:
  nops(S);
               3

> map(length,S);

              [951897, 1902263, 1902263]

It should be doable something like this:

> f:=(1/2)*exp(-RootOf(erf(_Z)*b*exp(exp((-x+a)/b))+
    b*exp(exp((-x+a)/b))-2*exp((-x+a)/b))^2-
    RootOf(erf(_Z)*sqrt(Pi)*d+sqrt(Pi)*d-sqrt(2)*
    exp(-(1/2)*(-y+c)^2/d^2))^2)/Pi;

> data:= [[1,2],[1.5,0.3],[1.1,1.4],[2.1,1],[1,1.3],[0.3,1]];
   with(Optimization):
   LL:= add(eval(ln(f), {x=p[1],y=p[2]}), p=data);
   NLPSolve(LL,a=0..1,b=1..2,c=0..1,d=1..2);

[-17.0044824339629983, [a = 0., b = 2., c = 0., d = 2.]]

For many parameter ranges, I get errors such as 

Error, (in Optimization:-NLPSolve) complex value encountered

Presumably, it is hitting cases where Maple finds a complex value for one of the RootOfs.

This forum is for Maple.  We don't do Mathematica here.  You might try wolfram.com.

Actually, it's the series for 0 <= |z| < 1.

In this case, for the first few terms of the series you can use

asympt(f(z), z);

Or for the complete series in symbolic form,

simplify(subs(t=1/z,convert(f(1/t),FormalPowerSeries,t)));

Rather more challenging would be to get the Laurent series valid for 1 < |z| < 2. 

  It would be a 3D surface, not a curve.  I can't tell from what you wrote whether what you want to do is interpolation: find a function z = f(x,y) of a given form to fit your data points (x[i],y[i],z[i]) exactly; or regression: find the function of a given form that best fits your data points. But Fit in the Statistics package should actually do both of these: it's interpolation if the number of independent parameters is equal to the number of data points, regression if it's less.  Caution: despite its apparent attractiveness, interpolation of more than a small number of points is likely to give spectacularly poor results.