I don't know if what I'm about to say is the cause of your problem, but even if it isn't, it seems likely to me that it will be the cause of a future problem.

I notice that you have made great use of the concatenation operator || to construct names (symbols in Maple-lingo). It is dangerous to use programmatically constructed symbols (constructed by any means, || or otherwise) explicitly in procedure code because All programatically constructed symbols are global. That includes symbols constructed by ||, cat, nprintf, parse, and convert(..., name). It doesn't matter if the "stem" of the symbol is declared local.

[Edit: I edited the above to highlight the distinction between names and symbols as those terms are used by Maple. Maple's symbols are defined at ?symbol (pretty much any sequence of characters, with enclosure in special quotes needed in some cases). Maple's names are a superset of symbols: names = symbols union indexed symbols.]

Here are two procedures for this. The first, SolveFunctionalRecurrence, computes the first N iterations of a functional recurrence, with N being an integer input. The second, PlotFunctionalRecurrence, simply plots all functions produced by the first.

SolveFunctionalRecurrence:= proc(
   #recurrence, solved for the highest-index function, with that highest index being just
   #a name:
   Rec::And(typefunc(name, indexed(name)), 1 &under nops) = algebraic,

   #parameter values (in form p::name= v::numeric) and initial functional conditions
   #in form (P::name[k::integer](t::name)= ...): 
   ICs::set({name, And(typefunc(name, indexed(integer)), 1 &under nops)}= algebraic),

   N::nonnegint #number of iterations to perform
)
local 
   n:= op(op(0, lhs(Rec))), j, k, F:= table(),
   ord:= -min(
      map(
         f-> op(op(0,f)), 
         indets(eval(rhs(Rec), n= 0), 'typefunc'(op([0,0], lhs(Rec))['integer'])) 
      )
   )
;
   for k to N do 
      F[k]:= eval['recurse'](Rec, ICs union {n= k, seq(F[k-j], j= 1..min(k-1,ord))})
   od;
   [entries(F, 'nolist', 'indexorder')]
end proc:

PlotFunctionalRecurrence:= proc(
   #output of SolveFunctionalRecurrence:
   F::list(And(typefunc(name, indexed(integer)), 1 &under nops) = algebraic),
   trange::name= range(realcons) #name and plot range of continuous variable
)
local k, n:= nops(F);
   plot(
      rhs~(F), trange, 
      'color'= [seq('COLOR'('HUE', (.85)*(k-1)/n), k= 1..n)],
      _rest
   )
end proc:   

And here's an example of it's usage on your functional recurrence. However, the results produced by your example initial function exp(-t) were too boring to be worth plotting, so I changed the initial function to the standard Normal distribution PDF, exp(-t^2/2)/sqrt(2*Pi).

F:= SolveFunctionalRecurrence(
   P[n](t) = (diff(P[n-1](t),t) + (lambda+mu)*P[n-1](t) - lambda*P[n-2](t))/mu,
   {P[-1](t)= 0, P[0](t)= exp(-t^2/2)/sqrt(2*Pi), lambda= 1, mu= 2},
   9 #number of iterations
);

PlotFunctionalRecurrence(F, t= -6..3, gridlines= false);

Note that my procedure requires that the recurrence be input in a form solved for P[n](t) with all functional terms on the right side having a lower index, and the number of initial functional conditions must match the order of the recurrence (even if that means that some must be specified as the 0 function, as is the case here).

The following one-line procedure removes n randomly selected elements from list L:

RemoveRand:= (L::list, n::nonnegint)-> (N-> L[[combinat:-randcomb(N, N-n)[]]])(nops(L)); 

There may be other ways that'll be overall more efficient if, say, L has thousands of elements or the operation is repeated hundreds of times.

Any plot can be animated, but you need to specify a parameter that varies with time. You have three potential parameters (r, theta, U), but you haven't quite specified which is supposed to vary with time. Does this code give what you want? I made U the parameter that varied with time.

plots:-animate(
   plots:-implicitplot,
   [
      U = r*cos(theta), r= 0..2, theta= 0..2*Pi, 
      coords= polar, axiscoordinates= polar
   ],
   U= -2..2
);

@yonasmibi If you want to use indexing, so that all previously computed values are stored in a table (the default) or a Vector (sometimes more efficient than a table) then you should do something like

for i to n do
   k[i+1]:= k[i]*(-b*r[i]+a+1);
   r[i+1]:= r[i]*(b*e*k[i]-c+1)
od;

This is the typical pattern for evaluating recurrence formulas: The index, i.e. the i or i+1 in this case, should be larger on the left side of := than on the right. To do this with Vectors, they should be declared with n+1 elements:

k:= Vector(n+1);  r:= Vector(n+1);

Tables never have a limit on their number of entries. They can be declared as k:= table(), etc., but it's often not necessary.

If you only want to get the final values after n iterations, then it's not necessary to use any indexing or container storage. All that's needed is

to n do (K,R):= (K*(-b*R+a+1), R*(b*e*K-c+1)) od:

It's necessary to use multiple assignment in the above. If you tried consecutive assignment, then the new value of K would erroneously get used in the formula for R.

For the distribution that you present, it is trivial and common to generate a sample simply by using rand; there's no need to refer to a PDF.

Sampler:= k-> n-> (r-> ['r()'$n])(rand(0..k)):
d4:= Sampler(4):
d4(9);

In addition to what vv said (which is totally correct):

The saving to a file and using in another worksheet is totally a red herring (English idiom for "misleading clue"). You'd have the same trouble if you tried to access psi in the first worksheet.

When you call newprocedure, regardless of which worksheet you call it from, you'll need to assign its results:

(H, psi):= newprocedure(2)

It is totally irrelevant whether you use the same names H and psi or different names in the line above. In other words, the names that you use for the matrices outside the procedure have no connection to the names used inside the procedure.

Yes, there is now a command that can do that: Explore. But I don't think it exists in Maple 12.

You can change the vertex labels with GraphTheory:-RelabelVertices. The help page ?GraphTheory,GraphDrawingOptions describes how to change the color, border, and font of the vertices. You can change the plot positions of the vertices with GraphTheory:-SetVertexPositions. But I don't think that there's any direct way to print the labels outside the vertex. However, it should be fairly easy to do it indirectly by using GraphTheory:-GetVertexPositions combined with plots:-textplot, GraphTheory:-DrawGraph(..., showlabels= false), and plots:-display.

If your expression containing theta is expr, then do 

diff(subs(theta= theta(t), expr), t);

In other words, theta has been replaced by theta(t), which indicates that theta is a function of t.

Here is a simple closed-form answer that only uses real numbers:

U := n->
   exp(.238022513478555*n)
    * (1.70226106278932*sin(1.91811758317245*n)
       + 1.28091347285205*cos(1.91811758317245*n))
    + 6.71908652714795*exp(.622567261710995*n);

To get exact integer values, apply round.

I am working on an elementary expository derivation of the above formula for you. 

Your initial conditions are wrong. They should be U(-2) = 2, U(-1) =2, U(0) = 2.

What I present here is not the standard straightforward explanation; rather, I think it's the most-compact elementary explanation. Nonetheless, you should be able to find it in any good textbook for a full-year introductory single-variable calculus course.

We'll suppose that a(x) > 0 for all x. This is only done pedagogically to avoid imaginary numbers; it's not mathematically necessary.

Let y(x) = a(x)^b(x)  =>   ln(y(x)) = b(x)*ln(a(x)) . Differentiate both sides wrt x:

y'(x)/y(x) = b'(x)*ln(a(x)) + b(x)*a'(x)/a(x)  =>  y'(x) = y(x)*(b'(x)*ln(a(x)) + b(x)*a'(x)/a(x))  

=>  y'(x) = a(x)^b(x)*(b'(x)*ln(a(x)) + b(x)*a'(x)/a(x)). Now, rewrite "prime" notation as "diff" notation and you'll have what you posted.

A list of lists, all having the same number of elements, is called a listlist in Maple-lingo. It's a primitive form of matrix, and it's often a useful substitute for an official Matrix. If A is a listlist, then its transpose is

`[]`~(A[])

Thus, what you want to do can be done by

A:= [[0, 1, 2, 3, 4, 5], [0, 2, 4, 6, 8, 10]]:
writedata(test, `[]`~(A[]), integer);
0	0
1	2
2	4
3	6
4	8
5	10

I don't build web pages, so I've never done this, but I suspect that you could embed a Maple Player  , which is a freely redistributable app whose main purpose is allowing those who don't own Maple to view your worksheets.

Maple can solve nonlinear ODE BVPs both symbolically and numerically. It can solve linear ODE IVPs by Laplace transforms, but not nonlinear ones.

Please post an example BVP of the type that you work with and state whether you want a symbolic or numeric solution.