coupled equations with finite difference method

You need to think of your coupled equations as one equation.

Be n=$left(begin{array}{c}n_0n_1n_2end{array}right)$, you can write it as $$frac{d}{dt}n= A n$$ with $$A=left(begin{array}{ccc}-W_{01}(t)&K_{10} &0W_{01}(t)&-W_{12}(t) - K_{10}&K_{21}&W_{12}(t)&-K_{21}end{array}right).$$

For the constraint you could either derive a new set of equations by explicitly solving the constraint for $n_2=1-n_0-n_1$ and plug it in or if you want to be flexible with linear constraints of such type, you could do the following.

Linear Constraints

We can rewrite the linear constraint equation as $$Bn=l$$ with $$B=left(begin{array}{ccc}1&1&1end{array}right)$$ and $$l=1.$$

Now be $I_B$ a matrix of a basis of the image space of $B$ and $N_B$ a matrix of a basis of the null space of B, we can do a change of basis by writing $n$ as $$n=I_B n_{I}+ N_B n_{N}.$$ The matrices $I_B$ and $N_B$ can be obtained by using a singular value decomposition on $B$. Our new variables are $n_I$ and $n_N$.

Using this basis in the constraint equation leads to $$B I_B n_I+B N_B n_N=l.$$ By construction $B N_B n_N= 0$, we solve the linear system for $n_I$ $$B I_B n_I = l.$$

The same basis can be used in the ODE system $$frac{d}{dt}left(I_B n_I + N_B n_N right)= A left(I_B n_I + N_B n_N right)$$ $frac{d}{dt}I_B n_I=0$ because it is a constant of motion due to the constraint equation.

Now multiplying both sides with $N_B^T$ leads to the reduced final system $$frac{d}{dt} n_N = N_B^T A I_B n_I + N_B^T A N_B n_N .$$

You can apply your favorite time integrator on the vectorized version of the equation for the variables $n_N$, since we haven taken care of the constraints. Remember, the first term on the rhs $N_B^T A I_B n_I$ is a known function of $t$, since we solved already for $n_I$.

Once you have obtained a solution for $n_N$, you can use our change of basis to obtain the solution in $n_0$ to $n_2$.

Time Integration

The simplest method is for sure the Euler method, which is in this case leads to $$n_N^{i+1}=n_N^i+Delta t left(N_B^T A(t_i) I_B n_I + N_B^T A(t_i) N_B n_N^iright).$$

Due to the dependency of the next step from the previous one, you cannot simply build up a system and solve it. Additionally, since $W_{01}(t)$ is a function of time, it is possible to a have a stiff or non stiff ODE system. So a recommendation is difficult, a good overview of possible solvers and there application can be found here.

Also you mention too many initial conditions. You only need at one point in time, e.g. the initial conditions at $t=0$, $$n_N(0)=N_B^Tleft(left(begin{array}{c}n_0(0)n_1(0)n_2(0)end{array}right)-I_B n_Iright).$$

I don't know scilab, but it seems that they have time integrators for ODEs, too. So check their documentation how to use them.