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}$$