Given two machines, both of which have an exponential lifetime with mean $\frac{1}{2}$. There is a single repairman that can service machines at an exponential rate $1$. Assuming that the system starts in a state in which both machines are working, using a numerical method, estimate the probability that both machines are working at time 1 using forward Kolmogorov equations.
My attempt: Let $X(t) =$ # of machines that failed by time $t$. So we have 3 states: $0$, $1$ and $2$ machines failed. Order the states as $0$, $1$ and $2$, then the rate transition matrix Q is: $$Q =\pmatrix{-4 & 4 & 0\\ 1 & -3 & 2\\ 0 & 1 & -1}$$
Now, using the forward Kolmogorov equation, we would get one of the $9$ equations is:
$P'_{00}(t) = -4P_{00}(t) + P_{01}(t)$
Thus, we only need to estimate $P_{00}(1)$ from numerically solving the differential equation above. Using Euler's method with stepsize $h=1$, since the system is currently in state $0$, we have: $P_{00}(0) = 1$ and $P_{01}(0) = 0$. Thus, by Euler's method, $P_{00}(1)\approx P_{00}(0) + h(-4P_{00}(0) + P_{01}(0)) = 1 + 1(-4+0) = -3 < 0$.
My question: Why did I get the negative probability? Is it because the stepsize $h=1$ is too large? But if I set $h<1$, why does it work? I don't see the difference except the increase in the accuracy of our estimated probability.