Hi, try this: load(solve_rec)$ solve_rec(x[n+2]-s*x[n+1]+p*x[n],x[n],x[0]=2,x[1]=s)$ define(x[n],rhs(%))$ [x[0],x[1],x[2]]$ ratsimp(%); should result: [2,s,s^2-2*p] Regards Marco Verpelli