diff --git a/sources/globals.f90 b/sources/globals.f90
index 55c82bf..70fe04a 100644
--- a/sources/globals.f90
+++ b/sources/globals.f90
@@ -16,7 +16,7 @@ module globals
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
 
 
-  CHARACTER(10), parameter :: version='v0.18.00' ! version number
+  CHARACTER(10), parameter :: version='v0.18.01' ! version number
 
 
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
@@ -238,6 +238,7 @@ module globals
   INTEGER              :: pp_nsteps      =  1
   INTEGER              :: pp_nfp         =  1
   REAL                 :: pp_xtol        =  1.000D-06
+  INTEGER              :: axis_npoints   = 360
                                                          
   namelist / focusin / &
   IsQuiet       ,&
@@ -337,7 +338,7 @@ module globals
   Nmax          ,&
   sdelta        ,&
   case_optimize ,&
-  fdiff_delta        ,&
+  fdiff_delta   ,&
   exit_tol      ,&
   DF_maxiter    ,&
   DF_xtol       ,&
@@ -376,7 +377,8 @@ module globals
   pp_maxiter    ,&
   pp_nsteps     ,&
   pp_nfp        ,&
-  pp_xtol                            
+  pp_xtol       ,&
+  axis_npoints                      
 
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
   
@@ -538,7 +540,7 @@ module globals
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
   
   ! fieldline tracing
-  REAL, ALLOCATABLE    :: XYZB(:,:,:), ppr(:,:), ppz(:,:), iota(:)
+  REAL, ALLOCATABLE    :: XYZB(:,:,:), ppr(:,:), ppz(:,:), iota(:), axis_phi(:), axis_r(:), axis_z(:)
   INTEGER              :: tor_num, total_num, booz_mpol, booz_ntor, booz_mn
   LOGICAL              :: lboozmn = .false.
   INTEGER, ALLOCATABLE :: bmim(:), bmin(:)
diff --git a/sources/poinplot.f90 b/sources/poinplot.f90
index 224206a..800d734 100644
--- a/sources/poinplot.f90
+++ b/sources/poinplot.f90
@@ -145,7 +145,8 @@ END SUBROUTINE poinplot
 !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
 
 SUBROUTINE find_axis(RZ, MAXFEV, XTOL)
-  USE globals, only : dp, myid, ounit, zero, pp_phi
+  USE globals, only : dp, myid, ounit, zero, pi2, sqrtmachprec, &
+                      pp_phi, pp_xtol, axis_phi, axis_r, axis_z, axis_npoints
   USE mpi
   IMPLICIT NONE
 
@@ -155,11 +156,13 @@ SUBROUTINE find_axis(RZ, MAXFEV, XTOL)
 
   INTEGER, parameter   :: n=2
   INTEGER              :: ml,mu,mode,nprint,info,nfev,ldfjac,lr
+  INTEGER              :: iwork(5), ierr, astat, ifail, i
   REAL     :: epsfcn,factor
   REAL     :: fvec(n),diag(n),qtf(n),wa1(n),wa2(n),wa3(n),wa4(n)
+  REAL     :: rz_end(n), phi_init, phi_stop, relerr, abserr, work(100+21*N)
   REAL, allocatable :: fjac(:,:),r(:)
-  external :: axis_fcn
- 
+  external :: axis_fcn, BRpZ
+
   LR = N*(N+1)/2
   LDFJAC = N
   ml = n-1
@@ -195,6 +198,22 @@ SUBROUTINE find_axis(RZ, MAXFEV, XTOL)
         write(ounit,'("findaxis: info="I2", something wrong with the axis finding subroutine.")') info
      end select
   endif
+ 
+  ! save axis data
+  SALLOCATE(axis_phi, (1:axis_npoints), zero)
+  SALLOCATE(axis_r  , (1:axis_npoints), zero)
+  SALLOCATE(axis_z  , (1:axis_npoints), zero)
+  relerr = pp_xtol
+  abserr = sqrtmachprec
+  rz_end = RZ
+  do i = 1, axis_npoints
+     axis_phi(i) = pp_phi + (i - 1) * pi2 / axis_npoints
+     axis_r(i) = rz_end(1) ; axis_z(i) = rz_end(2)
+     ifail = 1
+     phi_init = axis_phi(i)
+     phi_stop = phi_init + pi2 / axis_npoints
+     call ode( BRpZ, n, rz_end, phi_init, phi_stop, relerr, abserr, ifail, work, iwork )
+  enddo
 
   return
 
diff --git a/sources/saving.f90 b/sources/saving.f90
index 34307e5..3c4d895 100644
--- a/sources/saving.f90
+++ b/sources/saving.f90
@@ -157,6 +157,7 @@ subroutine saving
   HWRITEIV( 1                ,   pp_ns         ,   pp_ns                         )
   HWRITEIV( 1                ,   pp_maxiter    ,   pp_maxiter                    )
   HWRITERV( 1                ,   pp_xtol       ,   pp_xtol                       )
+  HWRITEIV( 1                ,   axis_npoints  ,   axis_npoints                  )
 
   HWRITEIV( 1                ,   Nfp           ,   surf(plasma)%Nfp                     )
   HWRITERV( 1                ,   surf_vol      ,   surf(plasma)%vol                     )
@@ -346,6 +347,9 @@ subroutine saving
      HWRITERA( pp_ns, pp_maxiter+1,   ppr         ,  ppr(1:pp_ns, 0:pp_maxiter) )
      HWRITERA( pp_ns, pp_maxiter+1,   ppz         ,  ppz(1:pp_ns, 0:pp_maxiter) )
      HWRITERV( pp_ns              ,   iota        ,  iota(1:pp_ns)              )
+     HWRITERV( axis_npoints       ,   axis_phi    ,  axis_phi(1:axis_npoints)   )
+     HWRITERV( axis_npoints       ,   axis_r      ,  axis_r(1:axis_npoints)     )
+     HWRITERV( axis_npoints       ,   axis_z      ,  axis_z(1:axis_npoints)     )
   endif
 
   if (allocated(XYZB)) then