  We apply a computationally efficient technique to validate the global
structure of the pulsar magnetosphere. In this first of a series of studies, a
3D, computationally intensive, implicit Crank-Nicolson finite-difference scheme
is developed. The region of magnetic influence is evolved under the
approximation of force-free electrodynamics. The main objective of this paper
is to present our code and use it to demonstrate and verify the now widely -
accepted global features of a pulsar magnetosphere. Our results qualitatively
agree with previously developed time-dependent models for an oblique rotator.
In line with earlier studies, we also demonstrate that our simulations can run
for many stellar rotations. Once we extend our code, we believe that our
implicit approach can be extremely useful to investigate magnetospheres filled
with resistive plasma, develop better resolution current sheets and investigate
small scale microphysics of pair creation using particle-in-cell techniques.
