The IVP should be sufficient to find the solution at any time. But with operator or matrix equations it does get confusing. A good approach is to solve the equation(s) first for regular functions and plug them in, seeing if they work for operators as well.
in this derivation: [solution1], the author assumes that $[S_i,S_j]=i \hbar \epsilon_{ijk}S_k$ doesn’t only hold for $t=0$, but for all times, why would that be?
I will show you here: We know that $U(t)U^\dagger (t) = U^\dagger(t) U(t) = 1$. And it is the case that we can write $S_i(t) = U^\dagger(t) S_i(0) U(t)$ - you can see this directly when trying to go from the schrödinger to the <PERSON> picture, and also I will justify it in my answer to your second question below. So, assuming that at t=0 the commutation holds, we multiply by a bunch of $U^\dagger$ from the left, by $U$ from the right, and insert the identity in the middle in order to get the time evolved version of all three matrices:
$$S_i S_j - S_j S_i = i \hbar \epsilon_{ijk} S_k$$ $$\implies U^\dagger (t) S_i U(t)\,\, U^\dagger(t) S_j U(t) - U^\dagger (t) S_j U(t)\,\, U^\dagger(t) S_i U(t) = i \hbar \epsilon_{ijk} U^\dagger(t)S_k U(t)$$ $$\implies S_i(t) S_j(t) - S_j(t) S_i(t) = i \hbar \epsilon_{ijk} S_k(t)$$
In this solution, discussed on this site: [solution2][2], the OP assumes that $S_j(t)$ is of the form $\exp(-iHt/\hbar)S_j(0) \exp(iHt/\hbar)$ (and therefore actually already postulates the answer), why can we assume that?
This is always true of any operator satisfying the heisenberg equation of motion, as long as it doesn't explicitly depend on $t$:
$$\frac{d}{dt}A(t)=\frac{i}{\hbar}[H,A(t)]$$
$$\implies A(t) = e^{i t H/\hbar}A(0)e^{-i t H/\hbar}$$
You can check this by plugging it in. 526
And this also means that $A(t)=U^\dagger (t)A(0)U(t)$, which is just the form of the equation above. It works with any observable operator, and any Hamiltonian. Pretty amazing that you can just write out the solution like that. But the caveat is that this only goes so far... there are much more useful ways to write the solution, as you saw in your post where we can write it explicitly in terms of pauli matrices for example. So some might not consider it actually "solved" if it is just written in this form.
Similarly, the solution in the <PERSON> picture for any wavefunction can be written
$$|\psi (t) \rangle = e^{-itH/\hbar}|\psi(0) \rangle$$ 526