lambda <- 1.71;
mu <- 2.97;
delta_t <- 2;
L <- 10;

#preparations, fucker!
data <- c(-lambda,             mu,                  0,                  0,       0,
           lambda, -(lambda + mu),             2 * mu,                  0,       0,
                0,         lambda, -(lambda + 2 * mu),             3 * mu,       0,
                0,              0,                  0,             lambda, -4 * mu,
                1,              1,                  1,                  1,       1);

left_side <- t(matrix(t(data), 5, 5));
right_side <- c(0, 0, 0, 0, 1);

#linear algebra, fucker!
p_fin <- solve(left_side, right_side);

#initial conditions and preparations, fucker!
times <- seq(0, (L - 1) * delta_t, by = delta_t);
params <- c(lambda, mu);
initial_cond <- c(1, 0, 0, 0, 0);

#differential equations, fucker!
System <- function(times, initial_cond, params) {
  with (as.list(c(initial_cond, params)), {
    dp0 <- -params[1] * initial_cond[1] + params[2] * initial_cond[2]
    dp1 <-  params[1] * initial_cond[1] - (params[1] + 1 * params[2]) * initial_cond[2] + 2 * params[2] * initial_cond[3]
    dp2 <-  params[1] * initial_cond[2] - (params[1] + 2 * params[2]) * initial_cond[3] + 3 * params[2] * initial_cond[4]
    dp3 <-  params[1] * initial_cond[3] - (params[1] + 3 * params[2]) * initial_cond[4] + 3 * params[2] * initial_cond[5]
    dp4 <-  params[1] * initial_cond[4] - 4 * params[2] * initial_cond[5]
    list (c(dp0, dp1, dp2, dp3, dp4))
  });
}
library("deSolve", lib.loc="C:/Dev/R-3.1.1/library");
out <- ode(y = initial_cond, times = times, func = System, parms = params, method = "rk4");