1/* Filename improved.mac
2
3   ***************************************************************
4   *							         *
5   *                     <package name>                          *
6   *                <functionality description>                  *
7   *                                                             *
8   *          from: Computer Algebra in Applied Math.            *
9   *                   by Rand (Pitman,1984)                     *
10   *                Programmed by Richard Rand                   *
11   *      These files are released to the public domain          *
12   *            						 *
13   * Reverse engineered from file newimprv.bk1                   *
14   * David Billinghurst - Feb 2004                               *
15   *                                                             *
16   ***************************************************************
17*/ /*   This program uses recursive functions to find
18the transition curves in Mathieu's equation.  To call it,
19type:
20                        TC()
21*/
22
23tc():=(input(),sign:1,find(),if n > 0 then (sign:-1,find()))$
24input():=(n:read("enter transition curve number n"),
25      m:read("enter degree of truncation"))$
26find():=(remarray(b,e),delta:n^2/4,for i thru m do delta:delta+d(i)*e^i,
27     print("delta=",delta),print(" "))$
28a(j,k):=
29  if j < 0 or k < 0 then
30    0
31  else if j = 0 and k = n then
32    1
33  else if j = 0 then
34    0
35  else if k = n then
36    0
37  else if numberp(b[j,k]) then
38    b[j,k]
39  else if k = 0 then
40    b[j,k]:(-a(j-1,2)/2-sum(d(i)*a(j-i,0),i,1,j))/(n^2/4)
41  else
42    b[j,k]:(-(a(j-1,k-2)+a(j-1,k+2)+sign*a(j-1,2-k))/2
43      -sum(d(i)*a(j-i,k),i,1,j))/((n^2-k^2)/4)
44$
45
46d(j):=
47  if numberp(e[j]) then
48    e[j]
49  else if n = 0 then
50    e[j]:-a(j-1,2)/2
51  else
52    e[j]:-(a(j-1,n-2)+a(j-1,n+2)+sign*a(j-1,2-n))/2
53$
54