All you need to do is solve your recursion and then sample using the resulting probabilities.

So for example in R you might assume $p(0)=0$ and use the following

```
max_n <- 20 # number of probabilities to calculate
p <- rep(0, max_n)
rec_p <- function(x1, steps=20){
p[1] <- x1
p[2] <- p[1]^2
for (i in 3:steps){ p[i] <- p[i-1]^2 + p[i-2]^3 }
p
}
err_rec_p <- function(x){ abs(sum(rec_p(x, max_n)) - 1) }
best_p1 <- optimize( f=err_rec_p, interval=c(0,1), tol = 1e-9 )$minimum
probabilities <- rec_p(best_p1, max_n) / sum(rec_p(best_p1, max_n))
```

to get the probabilities

```
[1] 5.005662e-01 2.505666e-01 1.882088e-01 5.115401e-02 9.283564e-03
[6] 2.200409e-04 8.485179e-07 1.137392e-11 6.110475e-19 1.471774e-33
[11] 2.281524e-55 3.188038e-99 1.187613e-164 3.240190e-296 0.000000e+00
[16] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
```

and then sample using these probabilities

```
samplesize <- 75
sample(1:max_n, samplesize, replace=TRUE, prob=probabilities )
```

which might for example give the following values

```
[1] 1 1 2 3 1 3 4 2 2 1 1 1 2 1 3 1 2 5 1 3 3 1 2 1 1 1 1 1 3 1 1 2 1 1 3 2 3 1
[39] 2 1 3 2 3 2 2 3 1 1 2 2 1 3 1 1 1 1 1 2 2 1 3 1 1 1 2 1 1 3 1 3 1 3 1 1 1
```