It can easily be done using for instance the freely available program Pari/GP. Here is a short excerpt of an example session.
We define the linear recurrence $x_{k+1} = B x_k + A x_{k-1} \quad $
For an example we set $[A,B]=[1,4] \quad $, then produce testdata beginning with $[x_1,x_2]=[1,2] \quad $. This is the vector of example data $[9, 38, 161, 682, 2889, 12238] \quad $ for which we want to find the recurrence-coefficients A and B (and possibly more if we do not know the order of recurrence beforehand).
Given for instance 6 numbers which are consecutive iterations, we could restate the problem as a 3x3-matrix equation: $ \small \begin{array} {} & \cdot & \begin{bmatrix} A \\ B \\ C \end{bmatrix} \\ \begin{bmatrix} 9 & 38 & 161 \\ 38&161&682 \\ 161&682&2889 \end{bmatrix} & = & \begin{bmatrix} 682 \\ 2889 \\ 12238 \end{bmatrix}\end{array}$ If the lhs-matrix has smaller rank than 3, then we need only two coefficients $[A,B] \quad $ instead of three $[A,B,C] \quad $.
The Pari/GP code to solve the problem:
\\pari/gp \\ leading ? indicates a command to Pari/GP, \\ leading % and a number indicates the (numbered) answer \\ we introduce the given vector of numbers ? v= [9, 38, 161, 682, 2889, 12238] %141 = [9, 38, 161, 682, 2889, 12238] \\ of six consecutive values we can generate a 3x4-matrix ? M=matrix(3,4,r,c,v[r+c-1]) %142 = [9 38 161 682] [38 161 682 2889] [161 682 2889 12238] \\ the rank of the matrix gives the order of the linear recurrence ? rk=matrank(M) %143 = 2 \\ so we'll have two coefficients and need a 2x2 and a 2x1 matrix to solve ? MX=matrix(rk,rk,r,c,M[r,c]) %144 = [9 38] [38 161] ? MY=matrix(rk,1,r,c,M[r,rk+1]) %145 = [161] [682] \\ the following command solves the problem: ? AB = MX^-1 * MY %146 = [1] [4]
The resulting matrix contains the coefficients A and B.
As a complete function this could be written:
{solvelinrec(v)=local(rs,cs,M,rk,MX,MY,AB); rs=#v \ 2; cs = rs+1 ; M=matrix(rs,cs,r,c,v[r+c-1]); rk=matrank(M); MX=matrix(rk,rk,r,c,M[r,c]); MY=vectorv(rk,r,M[r,rk+1]); AB = matsolve(MX,MY); \\ MX^-1 * MY ; return(AB); }
and then
? solvelinrec([9, 38, 161, 682, 2889, 12238]) %173 = [1, 4]~
or, receiving a order-three-recurrent sequence:
? solvelinrec([4, 5, 13, 11, 46, 11, 184, -94, 841, -1033, 4303, -7594]) %178 = [1, 4, -1]~
thus in this case we have $x_k = -1 x_{k-1} + 4 x_{k-2} + 1 x_{k-3}$