We consider the poisson equation where the functions involved are periodic including the solution function. Let $R=[0,1]{ imes}[0,l]{ imes}[0,1]$ be the region of interest and let $phi$(x,y,z) be an arbitrary periodic function defined in the region R such that $phi$(x,y,z) satisfies $phi$(x+1, y, z)=$phi$(x, y+1, z)=$phi$(x, y, z+1)=$phi$(x,y,z) for all x,y,z. We describe a very simple method for solving the equation ${ abla}^2u(x, y, z)$ = $phi$(x, y, z) based on the cubic spline interpolation of u(x, y, z); using the requirement that each interval [0,1] is a multiple of the period in the corresponding coordinates, the Laplacian operator applied to the cubic spline interpolation of u(x, y, z) can be replaced by a square matrix. The solution can then be computed simply by multiplying $phi$(x, y, z) by the inverse of this matrix. A description on how the storage of nearly a Giga byte for $20{ imes}20{ imes}20$ nodes, equivalent to a $8000{ imes}8000$ matrix is handled by using the fuzzy rule table method and a description on how the shape preserving property of the Laplacian operator will be affected by this approximation are included.