This one, of course, is easy, and rsolve will give you the solution.

> rsolve({N(t+1) = (1+b)*(1-d)*N(t), N(0) = N0}, N(t));

N0*(1-d+b-b*d)^t

If you assigned numerical values to N0, b and d, you can compute and plot points as follows.

 > n[0]:= N0;
    for i from 1 to stop_time do n[i]:= (1+b)*(1-d)*n[i-1] end do:
    plots[pointplot]([seq]([i,n[i]], i=0..stop_time);
   

Caution: only slightly more complicated recursions with symbolic parameters, or even just rational numbers rather than floats, are likely to crash Maple when doing this sort of iteration.

I think we may need a lot more explanation before we can understand exactly what you're trying to do.  What is your equation?  What is alpha?  Where are you trying to  approximate the solution?

You might also try

> plot([signum(cos(t))*abs(cos(t))^3, signum(sin(t))*abs(sin(t))^3, t = 0 .. 2*Pi]);

I'm not sure exactly what Maple is doing in this case, but here's how I would do it by hand.  Note that
sum(binomial(k,i)*t^i, i=0..k) = (1+t)^k (this is the binomial theorem).  Call this f(t).  Differentiating term by term, t*f'(t) = sum(i*binomial(k,i)*t^i,i=0..k).  Evaluate this at t = 1/(k-1) to get your result.

A truly one-dimensional plot would be rather strange, so you imbed the line in the plane, usually as the x axis.  Thus you might try:

> plots[pointplot]([seq]([x,0], x = [-2.5, -2,.5,1.5]), symbolsize=20,symbol=solidcircle,colour=red);

The "symbolic" version of pdsolve does not accept boundary conditions.  For the numeric version, these boundary conditions would be written as

 {u(x,0) = f(x), D[1](u)(x,0) = g(x)}

See the help page ?pdsolve,numeric.

Compute the values a[2,j] for whatever j values you're going to use, and then

plots[pointplot]([seq([1/n, a[2,2*n]/a[2,2*(n-1)]], n = whatever .. whatever)]);

Somewhat related questions were discussed recently here: www.mapleprimes.com/forum/howperformimmediateoutputprocedure, but AFAIK the print command always does immediate output.  What, for example, do you get with this?

> f:= proc()
local i, ti;
for i from 0 to 10 do
  print(2*i,"seconds");
  ti:= time();
  while time()<ti+2 do od:
od:
end proc;
f();

If you want a one-liner,

> OddEven:= L -> seq(L[i],i=map(op,[selectremove](type,[$1..nops(L)],odd)));

Warning: if you submit this, your prof is likely to guess that it's not your own work.

> map(`*`, y1, r-1);

If those are your differential equation and boundary conditions, the solution is simply F(y) = 1.  But I think you want F(-1) = -1.  Here's how I might proceed.
 

Your differential equation:

> de:= diff(F(s,R),s$4)-R*(diff(F(s,R),s)*diff(F(s,R),s$2)-F(s,R)*diff(F(s,R),s$3));   

Make it into a series in powers of R

> series(eval(de, F(s,R)=add(f[k](s)*R^k,k=0..5)),R);

Extract individual differential equations for each power up to R^5.

> des:= [seq(coeff(%,R,k),k=0..5)];

Define a function to make boundary conditions for f[k].

> bc:= k -> (f[k](-1)=`if`(k=0,-1,0),f[k](1)=`if`(k=0,1,0),D(f[k])(-1)=0,D(f[k])(1)=0);

Solve the boundary value problems for each power.

> for k from 1 to 6 do
   f[k-1]:= unapply(rhs(dsolve({des[k],bc(k-1)})),s)
 end do;

Here is a plot of the solutions for R=0 (in red) and R=20 (in blue).

> plot([f[0](s), add(f[k](s)*20^k, k=0..5)], s = -1 .. 1, colour=[red, blue]);

You seem to be insisting on doing things the hard way.  Moreover, you're using the old and deprecated linalg package rather than LinearAlgebra.  But anyway:
 

> PL := convert(P, list);
   vars:= map(parse, [$"a".."p"]);

Get some arbitrary values for v and _t[1] to _t[7], and get the variable values:

> vals:= subs(map(t -> (t = rand(0..100)()), indets(PL)),zip(`=`,vars,PL));

> subs(vals, Matrix(4,4,vars));

The first thing you have to do is tell us what you're trying to do.  What are k₁, k₂, ..., k_n supposed to be?  I can't make any sense of [1, 4, 11, 31, 83, 227, 616, 1674], and neither can the Encyclopedia of Integer Sequences or Maple's guessgf. 
 

Do you mean this?

> series(S, x);

-1/x+1-(1/4)*x+(1/128)*x^3-(21/1024)*x^4-(109/12288)*x^5+O(x^6)

 Or perhaps this?

> convert(S, FormalPowerSeries, x);

-1/x+Sum(Sum(2*x^k/alpha^(k+1), _alpha = RootOf(3072-1536*_Z+576*_Z^2-160*_Z^3+35*_Z^4)), k = 0 .. infinity)

Certainly not with evalf: that is for floating-point evaluation.  I would have thought that sum would do it, or else something from the SumTools package.  But none of them seems to work in this case.