t = var('t')
(k_0, k_1, k_2) = var('k_0, k_1, k_2, s_0, e_0')
c_p = function('c_p', t)
c_s = function('c_s', t)
c_e = function('c_e', t)
c_se = function('c_se', t)
eq_0 = -diff(c_p, t) +  +k_2*c_se
eq_1 = -diff(c_s, t) +  -k_0*c_s*c_e +k_1*c_se
eq_2 = -diff(c_e, t) +  -k_0*c_s*c_e +k_1*c_se +k_2*c_se
eq_3 = -diff(c_se, t) +  +k_0*c_s*c_e -k_1*c_se -k_2*c_se
eq_4 = c_e(0) == e_0
eq_5 = c_s(0) == s_0
eq_6 = c_se(0) == 0
eq_7 = c_p(0) == 0
desolve_system([eq_0, eq_1, eq_2, eq_3, eq_4, eq_5, eq_6, eq_7], [c_p, c_s, c_e, c_se])