This article focuses on developing and examining several numerical algorithms used to construct higher order Taylor methods for approximating the solution of a system of first order initial value problems. Some numerical test cases to demonstrate the efficiency of these algorithms are presented. The numerical results have shown to be consistent with the exact results.