R - Picard Method for General Polynomial

Question

I am currently writing a program in R to find solutions of a general polynomial difference equation using Picard's method. For an insight in the mathematics behind it (as math mode isn't available here): https://math.stackexchange.com/questions/2064669/picard-iterations-for-general-polynomials/2064732

Now since then I've been trying to work with the Ryacas package for integration. However I ran into trouble trying to work with the combination of expression and integration function.

library(Ryacas)

degrees = 3
a = c(3,5,4,6)
x0 = -1
maxIterations(10)


iteration = vector('expression', length = maxIterations)
iteration[1] = x0

for(i in 2:maxIterations){

  for(i in 1:degrees){

  exp1 = expression( a[i] * iteration[i-1] ^ i) 

  }

iteration[i] = x0 + Integrate(exp1, t)
}

but this results in

"Error in paste("(", ..., ")") : cannot coerce type 'closure' to vector of type 'character'"

and exp1 = expression(a[j] * iteration[i-1]^j) instead of an actual expression as I tried to achieve. Is there anyway I can make sure R reads this as a real expression (i.e. for example 3 * ( x0 ) ^ j for i = 2)?

Thanks in advance!

Edit: I also found the Subst() function, and currently trying to see if anything is fixable using it. Now I am mainly struggling to actually set up an expression for m coefficients of a, as I can't find a way to create e.g. a for loop in the expression() command.