Exploitation of geothermal resource results the decrease of fluid pressure in geothermal reservoir. In production process, the analysis of reservoir condition is made by observing the pressure state. When the fluid is pumped into the reservoir, the value of pressure varies with time which is called as transient state. The plot of pressure versus time will create a curve with a certain slope. From the graph of the pressure the reservoir condition can be analyzed. The solution of the pressure transient is identified as exponential integral equation Ei(x). When the input to the function is really small for example at x < 0.01, the equation will form an asymptotic curve. Analytical solution involves logarithm natural and Euler constant (γ). In this paper we try to approach the solution of exponential integral equation by numerical integration. The objective of this study is to make a numerical model of the pressure change in a geothermal reservoir and to compare the result between numerical method and analytic. There are two methods used in this study, first is Picard- McLaurin iteration to solve the ordinary differential equations (ODE), and the second is trapezoidal integration to calculate the function of Ei(x). The modeling shows that the result of the calculation with the numerical method matched with the analytic with the range of error between 0.0008 to 4.5 % for drawdown test and 0.19 to 7.7 % for buildup test.