I tried the following code...... un := [0, 0, 0, 0]: un[1] := tanh(k)*tanh(k*n): for i from 2 to 6 do uns[i-1] := eval(un[i-1], t = s): und[i-1] := diff(uns[i-1], s): unp[i-1] := eval(uns[i-1], n = n+1): unm[i-1] := eval(uns[i-1], n = n-1): un[i] := un[i-1]-(int((und[i-1]-(1-uns[i-1])^2)(unp[i-1]-unm[i-1]), s = 0 .. t)): end do: