This paper describes numerical methods and the corresponding maximum-norm error estimates for a chemotaxis model proposed by Kepler. Upon detecting pathogens, immune cells secrete soluble factors that attract other immune cells to the site of the infection. The motion of the model cells is stochastic, but biased toward the gradient of one or more of the soluble factors. The soluble factors are modeled by a system of reaction-diffusion equations with sources that are centered on the cells. Previously, I presented a first order splitting in time for solving the reaction-diffusion-stochastic system numerically. The diffusion, reaction, and stochastic differential equations can be approximated separately to first order in the supremum norm. The three-dimensional domain is discretized using finite elements, and the diffusion is solved using a backward Euler scheme combined with multigrid. The reaction is solved using a simple semi-implicit first order scheme. The stochastic differential equations are given by a Langevin process which can be simulated exactly. The paper concludes by demonstrating first order convergence of the entire simulation and providing a sample simulation of the reaction-diffusion-stochastic system.