What you're running into is called "premature evaluation".  Your function f works fine with a numerical input, but not with a symbolic input.  Now sum, like most Maple procedures, starts off by evaluating its inputs, so when you write

> sum(f(t), t = -2 .. 2);

the first thing that happens is that the function f is called with the symbolic input t, and this causes an error as you have seen.  There are three basic strategies to get around this, and other posters have suggested examples of all three:

1) Use add instead of sum, because add has special evaluation rules and does not evaluate its first input.

2) Delay evaluation using quotes.

3) Make a procedure that will work for symbolic inputs, e.g. using piecewise instead of if.

Along the lines of #3, I might suggest using type-checking to make a procedure that will return unevaluated if given non-numeric inputs.

> f := proc(i) 
      if not type(i, numeric) then 'procname'(i)
      elif i > 0 then i
      else -i
      end if
   end proc;

Using reduction of order, since y[1](x) = (1-x^2)^(n/2) is one solution, the general solution is G(x) (1-x^2)^(n/2) where g(x) = G'(x) satisfies a first order differential equation

> de := (1-x^2)*diff(y(x),`$`(x,2))-2*x*diff(y(x),x)+(n*(1+n)-n^2/(1-x^2))*y(x);

   eval(de,y(x)=(1-x^2)^(n/2)*G(x));

   normal(%);

-(1-x^2)^(1/2*n)*(diff(G(x),`$`(x,2))*x^2+2*n*x*diff(G(x),x)+2*x*diff(G(x),x)-diff(G(x),`$`(x,2)))

> de1:= eval(op(-1,%),diff(G(x),x)=g(x));

de1 := diff(g(x),x)*x^2+2*n*x*g(x)+2*x*g(x)-diff(g(x),x)

> dsolve(de1);

g(x) = (x-1)^(-n-1)*(x+1)^(-n-1)*_C1

 > convert(rhs(%),FormalPowerSeries,x) assuming n::posint;

Sum(-_C1*(-1)^n*pochhammer(1+n,k)/k!*x^(2*k),k = 0 .. infinity)

> G:= int(%,x);

G := Sum(-_C1*(-1)^n*pochhammer(1+n,k)/k!*x^(2*k)*x/(1+2*k),k = 0 .. infinity)

Since the series of (1-x^2)^(n/2) is

Sum(pochhammer(-1/2*n,k)/k!*x^(2*k),k = 0 .. infinity)

we get our second solution as a power series

sum(a[k]*x^(2*k+1),k=0..infinity)

(that should be: sum(a[k]*x^(2*k+1),k=0..infinity)) where a[k] is

Sum((-1)^(n+1)*pochhammer(1+n,j)/j!/(1+2*j)*pochhammer(-n/2,k-j)/(k-j)!,j=0..k)

 Maple evaluates this sum as

-(-1)^n*pochhammer(-1/2*n,k)/k!*hypergeom([1/2, -k, 1+n],[3/2, 1+1/2*n-k],1)

which seems to be OK if n is odd, but not if n is even. 

I don't understand what you mean by "split" a function.  Can you explain, and/or give an example?

You can use normal axes, but add the tickmarks manually, like this:

> with(plots):
 for n from 0 to 19 do
   ticks:= [seq(evalf(j-Pi/10*n),j=ceil(Pi/10*n) .. floor(Pi/10*n+2*Pi))];
   frame[n]:= display([plot(sin(x+Pi/10*n),x=0..2*Pi,xtickmarks=[]),
     seq(plot([[t,0],[t,-0.02]],colour=black),t=ticks),
     seq(textplot([ticks[j],-0.1,j-1+ceil(Pi/10*n)]),j=1..nops(ticks))]);
 end do:
 display([seq(frame[n],n=0..19)],insequence=true);

To see this move, I think you have to double-click on it twice.

I assume you mean

diff(q,t)=C*diff(q,z,z)

If C > 0, this is the heat equation with time reversed.  If you want closed-form solutions:

> pde:= diff(q(z,t),t) = C*diff(q(z,t),z,z);
   pdsolve(pde, build); 

q(z,t) = exp(_c[1]^(1/2)*z)*_C3*exp(C*_c[1]*t)*_C1+1/exp(_c[1]^(1/2)*z)*_C3*exp(C*_c[1]*t)*_C2

As for Fourier series: note that the closed-form separation-of-variables solution will involve sines
or cosines in z, for negative _c[1].  For example, if your boundary conditions are q(0,t) = q(L,t) = 0,
the appropriate solution involves sin(Pi*n*z/L) for integers n.

> eval(pde, q(z,t)= sin(Pi*n*z/L)*exp(a*t));

sin(Pi*n*z/L)*a*exp(a*t) = -C*sin(Pi*n*z/L)*Pi^2*n^2/L^2*exp(a*t)

> solve(identity(%,z),a);

-C*Pi^2*n^2/L^2

(that should be -C*Pi^2*n^2/L^2)

So the Fourier series solution is

sum(c[n]*sin(Pi*n*z/L)*exp(-C*Pi^2*n^2/L^2*t),n = 1 .. infinity)

If you have an initial condition q(z,0) = g(z), you can get the coefficients c[n] from the Fourier series for g(z).

The trouble is, for typical initial conditions the series is likely not to converge for positive t.

The exponential generating function for the left side is (for sufficiently small s)

> G := sum(sum(k^n*x^k*s^n/n!, n=0..infinity), k=1..infinity);

G := -1/(x*exp(s)-1)*x*exp(s)

That should be: G := -1/(x*exp(s)-1)*x*exp(s)

So you want to prove that the n'th derivative of G wrt s, evaluated at s=0, is the right side of your equation.  In fact, since G depends only on x*exp(s), this amounts to showing that

(1-exp(s))^(n+1)*Diff(1/(1-exp(s))*exp(s),`$`(s,n)) = Sum(Sum((-1)^(m+1)*binomial(n+1,m-1)*(r-m+1)^n,m = 1 .. r)*exp(r*s),r = 1 .. n)

Now if the left side is P[n](exp(s)) , we have P[1](exp(s)) = exp(s) (which does agree with the right side for n=1), and

P[n+1](exp(s)) = (1-exp(s))^(n+2)*Diff((1-exp(s))^(-n-1)*P[n](exp(s)),s)

>  (1-exp(s))^(n+2)*diff((1-exp(s))^(-n-1)*P[n](exp(s)),s);
   P[n+1](exp(s)) = map(factor,collect(simplify(%),[P[n],D(P[n])]));
P[n+1](exp(s)) = exp(s)*(n+1)*P[n](exp(s))-D(P[n])(exp(s))*exp(s)*(-1+exp(s))

Or with exp(s) = u:

> subs(exp(s)=u,%);

P[n+1](u) = u*(n+1)*P[n](u)-D(P[n])(u)*u*(-1+u)

By induction on n, it's easy to see that P[n] will be a polynomial of degree n with constant term 0.   So let's write

P[n](u) = sum(a[n,r]*u^r,r=1..n)

> map(factor,collect(eval(u*(n+1)*P[n](u)-D(P[n])(u)*u*(-1+u),P[n]=(u->a[n,r]*u^r)),u));

a[n,r]*u^r*(n+1-r)*u+a[n,r]*u^r*r

So a[n,r] must satisfy the recursion

 a[n+1,r]=(n-r+2)*a[n,r-1]+r*a[n,r]

That should be: a[n+1,r]=(n-r+2)*a[n,r-1]+r*a[n,r]

At this point I have some work to do, so I'll leave it to you to show that Sum((-1)^(m+1)*binomial(n+1,m-1)*(r-m+1)^n,m = 1 .. r) satisfies this recursion.

The property should apply to the function sin, not to the expression sin(x).

> is(sin, continuous);

true

> is(sin, monotonic);

false

I don't know any way to say in Maple that a function is monotonic on some interval, rather than everywhere.
Can we add the assumption that sin is monotonic?

> assume(sin, monotonic);

Error, (in assume) cannot assume on an assigned name
 

However, this works:

>  `property/object`[sin]:= AndProp(op(`property/object`[sin]), monotonic);
    forget(is);
    is(sin, monotonic);

true

Most serious Maple programmers do not edit their code in a Maple worksheet. 

Use your favourite editor, produce a text file, and then read it in to Maple with the read command.

> pde:= diff(C(r, t), t) = Di*(diff(r^2*(diff(C(r, t), r)), r))/r^2-lambda*C(r, t)+Beta;
  sol:= dsolve(eval(pde, C(r,t) = C(r)));
sol := C(r) = 1/r*sinh((lambda/Di)^(1/2)*r)*_C2+1/r*cosh((lambda/Di)^(1/2)*r)*_C1+Beta/lambda

Note that for a limit to exist as r -> 0, you need _C1 = 0 (this can be seen using limit, but it's simple enough to see by inspection).

> limit(rhs(sol), r=0, right);

signum(_C1)*infinity

> sol := eval(sol, _C1=0);
  limit(rhs(sol), r=0, right);

(_C2*lambda*(lambda/Di)^(1/2)+Beta)/lambda

> solve(eval(rhs(sol), r=a), _C2);
  sol := eval(sol, _C2=%);
  sol := C(r) = -1/r*sinh((lambda/Di)^(1/2)*r)*Beta*a/sinh((lambda/Di)^(1/2)*a)/lambda+Beta/lambda

conjugate(X)

Another way is to do a parametric plot.  For example:
 

> f := t -> Int(sin(x*t)/(1+x^2), x=0..Pi);
   plot([ f(t), t, t = 0 .. 1]);

I'd generally prefer not to use the Explicit or _EnvExplicit option, which can produce horrendously large expressions for a quartic equation (bad enough in this particular case, but I've seen much worse).

> solutions := solve({e1,e2,e3});

solutions := {a = 33/2*RootOf(2*_Z^2-4*_Z+1,label = _L1)-6, x = 3*RootOf(2*_Z^2-4*_Z+1,label = _L1), y = 3*RootOf(2*_Z^2-4*_Z+1,label = _L1)}, {a = 15/13*RootOf(81+_Z^4-14*_Z^3-78*_Z+65*_Z^2)^2+51/26-3/26*RootOf(81+_Z^4-14*_Z^3-78*_Z+65*_Z^2)^3-18/13*RootOf(81+_Z^4-14*_Z^3-78*_Z+65*_Z^2), x = 10/13*RootOf(81+_Z^4-14*_Z^3-78*_Z+65*_Z^2)^2-1/13*RootOf(81+_Z^4-14*_Z^3-78*_Z+65*_Z^2)^3-25/13*RootOf(81+_Z^4-14*_Z^3-78*_Z+65*_Z^2)+30/13, y = RootOf(81+_Z^4-14*_Z^3-78*_Z+65*_Z^2)}

Note that there are two solution families, one involving a quadratic and one involving a quartic.

> P1:= indets(solutions[1],RootOf):
   rts1:= fsolve(op([1,1],P1));

rts1 := .2928932188, 1.707106781

> map(r -> subs(op(P1)=r, solutions[1]), {rts1});

   {{a = -1.167261890, x = .8786796564, y = .8786796564}, {a = 22.16726189, x = 5.121320343, y = 5.121320343}}

> P2:= indets(solutions[2],RootOf):
  rts2:= fsolve(op([1,1],P2));

rts2 := ``

The quartic has no real roots.  If you were interested in complex solutions as well:

>  rts2:= fsolve(op([1,1],P2), _Z, complex);

rts2 := .5845240526-1.155216020 I, .5845240526+1.155216020 I, 6.415475947-2.676840665 I, 6.415475947+2.676840665 I

> map(r -> subs(op(P2)=r, solutions[2]), {rts2});

{{a = .2535721579+0.*I, x = .584524053-1.155216020*I, y = .5845240526+1.155216020*I}, {a = .2535721579+0.*I, x = .584524053+1.155216020*I, y = .5845240526-1.155216020*I}, {a = 17.74642784-.2e-8*I, x = 6.41547595+2.676840659*I, y = 6.415475947-2.676840665*I}, {a = 17.74642784+.2e-8*I, x = 6.41547595-2.676840659*I, y = 6.415475947+2.676840665*I}}

For a simple SI model:

> model:= {D(i)(t) = beta*i(t)*(n-i(t)), i(0) = i0};
   solution:= dsolve(model);

solution := i(t) = -i0*n/(-i0-exp(-t*n*beta)*n+exp(-t*n*beta)*i0)

 > plot(eval(rhs(solution), {beta = 0.02, i0 = 1, n = 100}), t=0..5);

Why not just solve directly for S?

> mm:=(101/100)*ln(10)-ln(S)+10-S = 5*t;
   S1:= solve(mm, S);
 

S1 := exp(-LambertW(exp(10-5*t+101/100*ln(10)))+10-5*t+101/100*ln(10))

And then for the difference between S1 and the solution from dsolve

> plots[odeplot](sol1, [t, S(t)-S1], t=0..8, colour=red);

Find the height of the peak, divide by two, find the points where the function has that value, and subtract.  Depending on the function, you may want to use fsolve rather than solve.