The inhomogeneous burnup equation is often used for describing the time evolution of nuclides' depletion in nuclear systems which have a significant nuclide migration effect. However, lots of burnup calculations codes only deal with the homogeneous cases instead of the inhomogeneous ones, among them there are a few codes that can work only when the inhomogeneous term of the equation is constant. Based on the condition that the inhomogeneous term can be approximated by finite-order Taylor expansion, two methods are introduced to solve the inhomogeneous burnup equation whose inhomogeneous term is time dependent. For the first method, the transmutation trajectory analysis method is used to decompose the connections between nuclides into linear chains, for one chain the analytical solution is derived strictly by using the Laplace transform. For the second method, a solution of the inhomogeneous equation in the form of summation of infinite matrix series is first derived, and then the sum function of the series is found. Furthermore, the different-order nearly-best rational approximation function of the sum function is found by using Carathéodory-Fejér method. The error between the sum function and the rational function fluctuates in a certain range without exceeding a limit value, while the maximum error decreases exponentially with the order of rational function increasing. By adopting the nearly-best rational approximation, the summation of infinite matrix series converts into a finite expansion of matrix fraction, which is much easier to deal with. These two methods are implemented in the burnup calculation code JBURN and numerical tests are done through using two examples. The first example is a small-scale matrix example and the result shows that the results from the two methods agree well in at least 6 decimal precision together with the results from the reference solution. The second example is a large-scale problem based on real nuclides' reaction database, and the result shows that less than 1% among all nuclides have a deviation larger than 10% between two methods, while about 8% nuclides have a deviation larger than 0.01% and the remaining ones have a deviation smaller than 0.01%. These results validate the correctness and accuracy for each of the two methods. Finally, this paper provides a possible implementation process for solving inhomogeneous burnup equations which have other time-dependent forms of inhomogeneous term.