Below is my code for numerical solving of PDE with Crank Nicolson scheme.

 ClearAll[L, M, \[CapitalDelta]t, \[CapitalDelta]x]
    L = Pi // N; M = 25;
    \[CapitalDelta]x = L/M;
    \[CapitalDelta]t = \[CapitalDelta]x^2; k = 1;
    U0[x_] = x (Pi - x);
    U[0, n_] := 0.;
    U[M, n_] := 0.;
    sol[0] = Table[U[j, 0] -> U0[j \[CapitalDelta]x], {j, 1, M - 1}];
    sol[n_] := 
     Module[{vars, eqns}, vars = Table[U[j, n], {j, 1, M - 1}]; 
      eqns = Table[
         U[j, n] - U[j, n - 1] == 
          1/2 ((k \[CapitalDelta]t )/\[CapitalDelta]x^2 (U[j + 1, n] - 
               2 U[j, n] + U[j - 1, n] + U[j + 1, n - 1] - 2 U[j, n - 1] +
                U[j - 1, n - 1])), {j, 1, M - 1}] /. sol[n - 1];
      Solve[eqns, vars][[1]]

The problem is, that I have to eveluate value of scheme for (Pi/4,2). That means i have to put 127 for n in code below. The problem is, that the program doesn't want to eveluate for such big numbers.

rez = Table[
  Table[{\[CapitalDelta]x j, \[CapitalDelta]t n, U[j, n]}, {j, 0, 
     M}] /. sol[n], {n, 127}]

But if you go bigger for n, like 100, 400, and more it gives:

$RecursionLimit::reclim: Recursion depth of 256 exceeded. >>
  • 1
    $\begingroup$ Then set $RecursionLimit to a bigger value, or even Infinity. You will have to be careful if you do $RecursionLimit = Infinity, though. $\endgroup$ May 18, 2013 at 17:14
  • 1
    $\begingroup$ @J.M. Actually, I think that setting $RecursionLimit = Infinity is never appropriate. I realize that I may be biased by my software development (rather than problem - solving) perspective, but it is all too easy to lose unsaved data by using this setting. $\endgroup$ May 18, 2013 at 17:22
  • $\begingroup$ @Leonid, of course, I usually do that as Block[{$RecursionLimit = Infinity}, (* stuff *)] instead of setting it globally. I did say that OP has to be careful; in particular, he should be sure that the recursion he has is guaranteed to terminate, but usually people are unable or unwilling to do such proofs. $\endgroup$ May 18, 2013 at 17:26
  • $\begingroup$ Thanks. What is the correct command for blocking? Block[{$RecursionLimit = Infinity}, rez] in my example? $\endgroup$
    – energyMax
    May 18, 2013 at 17:29
  • 2
    $\begingroup$ @J.M. It also depends on how recursion is used. If one uses it for algorithm which is deeply recursive (high recursion depth), then the right way to do this in Mathematica is to perform some form of tail call optimization. Note that tail call - optimized functions in Mathematica require modifications of $IterationLimit rather than $RecursionLimit, and the former are safe (it is safe to set $IterationLimit to Infinity). $\endgroup$ May 18, 2013 at 17:34

1 Answer 1


Ok, your case is the one which can be easily helped by memoization. Basically, it will serve to compute the set of values in sol in reverse order, de facto making it just a convenient device to turn your algorithm into an iterative one without changing it much.

All you have to do is to change your definition of sol as:

sol[n_] := sol[n] = Module[your-old-code]

and then, if you want say n=1200, first execute


This will automatically add definitions for specific values of sol for all n from 1 to 1200, and then you simply call



and you get this under a couple of seconds.

  • $\begingroup$ Works as charm, thanks! $\endgroup$
    – energyMax
    May 18, 2013 at 18:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.