diff --git a/.gitignore b/.gitignore
index 491f19f94ad20f3f9e7c4d6e0729a9c6fabfc3e7..0b32f91f21a82f38540d5e890279fb4e5a9c0278 100644
--- a/.gitignore
+++ b/.gitignore
@@ -19,6 +19,7 @@
 *.h5
 *.o
 *.a
+*.png
 logs/
 results/
 results_old/
diff --git a/matlab/create_gif.m b/matlab/create_gif.m
index 826ea6def012776a647ec57e97ccd545e6f2ed0b..f2445a7f9af68af4b4ce9442aa16ade29d80869a 100644
--- a/matlab/create_gif.m
+++ b/matlab/create_gif.m
@@ -13,7 +13,7 @@ else
 % Setup figure frame
 fig  = figure('Color','white','Position', [100, 100, 400, 400]);
     pcolor(X,Y,FIELD(:,:,1)); % to set up
-    colormap jet
+    colormap gray
     axis tight manual % this ensures that getframe() returns a consistent size
     if INTERP
         shading interp;
diff --git a/matlab/create_gif_3D.m b/matlab/create_gif_3D.m
index 188349dba19c3e7c75e625a5af0d940fa2a3ddc2..646080fc0920fddb4f0d48db3d81fa92ccb8cf2d 100644
--- a/matlab/create_gif_3D.m
+++ b/matlab/create_gif_3D.m
@@ -10,7 +10,6 @@ fig  = figure('Color','white','Position', [100, 100, 400, 400]);
     scale_z = max(abs(plt(Z(:,:,1))));
     plot3(plt(X)/scale_x,plt(Y)/scale_y,plt(Z)/scale_z,'.k','MarkerSize',MARKERSIZE);
     view(VIEW);
-    colormap jet
     axis tight manual % this ensures that getframe() returns a consistent size
     in      = 1;
     nbytes = fprintf(2,'frame %d/%d',in,numel(FRAMES));
diff --git a/matlab/default_plots_options.m b/matlab/default_plots_options.m
index ed53cc039f90ba647c6f2c417f7ad72e716e7578..56e500d37e083493255aef09c715982e40ac57ea 100644
--- a/matlab/default_plots_options.m
+++ b/matlab/default_plots_options.m
@@ -6,4 +6,5 @@ set(0,'defaultAxesFontSize',16)
 set(0, 'DefaultLineLineWidth', 2.0);
 line_colors = [[0, 0.4470, 0.7410];[0.8500, 0.3250, 0.0980];[0.9290, 0.6940, 0.1250];...
     [0.4940, 0.1840, 0.5560];[0.4660, 0.6740, 0.1880];[0.3010, 0.7450, 0.9330];[0.6350, 0.0780, 0.1840]];
-line_styles = {'-','-.','--',':'};
\ No newline at end of file
+line_styles = {'-','-.','--',':'};
+marker_styles = {'^','v','o','s'};
\ No newline at end of file
diff --git a/matlab/save_figure.m b/matlab/save_figure.m
index e71ee646951dc1664ee43b5f009d07f451298b36..c5a9c72c39868cd1b4f840b20075f971e1c77c58 100644
--- a/matlab/save_figure.m
+++ b/matlab/save_figure.m
@@ -1,5 +1,8 @@
 %% Auxiliary script to save figure using a dir (FIGDIR), name (FIGNAME) 
 %  and parameters
-FIGNAME = [BASIC.RESDIR, FIGNAME,FMT];
-saveas(fig,FIGNAME);
-disp(['Figure saved @ : ',FIGNAME])
\ No newline at end of file
+if ~exist([BASIC.RESDIR,'/fig'], 'dir')
+   mkdir([BASIC.RESDIR,'/fig'])
+end
+saveas(fig,[BASIC.RESDIR,'/fig/', FIGNAME,'.fig']);
+saveas(fig,[BASIC.RESDIR, FIGNAME,'.png']);
+disp(['Figure saved @ : ',[BASIC.RESDIR, FIGNAME,'.png']])
\ No newline at end of file
diff --git a/wk/analysis_2D.m b/wk/analysis_2D.m
index 1ffdb66d1705c94f5184869cc8f5dcfb5d9f0693..6c739bb4c5052b3fdae5a3523f85c80217f0bf98 100644
--- a/wk/analysis_2D.m
+++ b/wk/analysis_2D.m
@@ -1,4 +1,5 @@
 %% Load results
+% JOBNUM = 0; load_results;
 compile_results
 
 %% Retrieving max polynomial degree and sampling info
@@ -165,7 +166,7 @@ save_figure
 %%
 if 1
 %% Show frame in real space
-tf = 100; [~,it] = min(abs(Ts2D-tf)); [~,it5D] = min(abs(Ts5D-tf));
+tf = 1000; [~,it] = min(abs(Ts2D-tf)); [~,it5D] = min(abs(Ts5D-tf));
 fig = figure; FIGNAME = ['rz_frame',sprintf('_%.2d',JOBNUM)];
     subplot(221); plt = @(x) (((x)));
         pclr = pcolor((RR),(ZZ),plt(ne00(:,:,it))); set(pclr, 'edgecolor','none'); colorbar;
diff --git a/wk/fort.90 b/wk/fort.90
index c849db79e720ade275ed7a76d4d88da685d19064..f684ea3ed002a4a38469b7a7a93836c8a4e43276 100644
--- a/wk/fort.90
+++ b/wk/fort.90
@@ -1,7 +1,7 @@
 &BASIC
   nrun   = 100000000
   dt     = 0.01
-  tmax   = 200
+  tmax   = 300
   RESTART = .false.
 /
 &GRID
@@ -25,9 +25,9 @@
   write_phi     = .true.
   write_non_lin     = .true.
   write_doubleprecision = .true.
-  resfile0      = '../results/ZP/256x128_L_66_Pe_2_Je_1_Pi_2_Ji_1_nB_0.5_nN_1_nu_1e-02_DG_mu_5e-04/outputs'
-  rstfile0      = '../results/ZP/256x128_L_66_Pe_2_Je_1_Pi_2_Ji_1_nB_0.5_nN_1_nu_1e-02_DG_mu_5e-04/checkpoint'
-  job2load      = 1
+  resfile0      = '../results/ZP/256x128_L_66_Pe_2_Je_1_Pi_2_Ji_1_nB_0.4_nN_1_nu_1e-02_DG_mu_5e-04/outputs'
+  rstfile0      = '../results/ZP/256x128_L_66_Pe_2_Je_1_Pi_2_Ji_1_nB_0.4_nN_1_nu_1e-02_DG_mu_5e-04/checkpoint'
+  job2load      = 0
 /
 &MODEL_PAR
   ! Collisionality
@@ -44,7 +44,7 @@
   q_i     = 1
   eta_n   = 1
   eta_T   = 0
-  eta_B   = 0.5
+  eta_B   = 0.4
   lambdaD = 0
   kr0KH   = 0
   A0KH    = 0
diff --git a/wk/load_results.m b/wk/load_results.m
index 5fccf08ef6a800454e2452f411f99ca39804dca3..43bdbdeb1d508f8aef6373ca1528731c7b1e45f2 100644
--- a/wk/load_results.m
+++ b/wk/load_results.m
@@ -2,7 +2,8 @@
 filename = sprintf([BASIC.RESDIR,'outputs_%.2d.h5'],JOBNUM);
 disp(['Loading ',filename])
 % Loading from output file
-% CPUTIME   = h5readatt(filename,'/data/input','cpu_time');
+CPUTIME   = h5readatt(filename,'/data/input','cpu_time');
+DT_SIM    = h5readatt(filename,'/data/input','dt');
 if strcmp(OUTPUTS.write_moments,'.true.') 
 [Nipj, Pi, Ji, kr, kz, Ts5D, dt5D] = load_5D_data(filename, 'moments_i');
 [Nepj, Pe, Je                    ] = load_5D_data(filename, 'moments_e');
diff --git a/wk/parameters_ZP.m b/wk/parameters_ZP.m
index 4033fbd2410165c2b4c623a87e30df5dd6af52b9..4339223c3486054826a8097812b3a56a62e59bc3 100644
--- a/wk/parameters_ZP.m
+++ b/wk/parameters_ZP.m
@@ -1,13 +1,12 @@
 clear all;
 addpath(genpath('../matlab')) % ... add
-default_plots_options
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Set Up parameters
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% PHYSICAL PARAMETERS
 NU      = 1e-2;   % Collision frequency
 TAU     = 1.0;    % e/i temperature ratio
-ETAB    = 0.5;    % Magnetic gradient
+ETAB    = 0.4;    % Magnetic gradient
 ETAN    = 1.0;    % Density gradient
 ETAT    = 0.0;    % Temperature gradient
 MU      = 5e-4;   % Hyper diffusivity coefficient
@@ -20,12 +19,12 @@ JMAXE   = 1;     % Highest ''       Laguerre ''
 PMAXI   = 2;     % Highest ion      Hermite polynomial degree
 JMAXI   = 1;     % Highest ''       Laguerre ''
 %% TIME PARAMETERS 
-TMAX    = 200;  % Maximal time unit
+TMAX    = 300;  % Maximal time unit
 DT      = 1e-2;   % Time step
 SPS2D   = 2;      % Sampling per time unit for 2D arrays
 SPS5D   = 0.1;    % Sampling per time unit for 5D arrays
 RESTART = 0;      % To restart from last checkpoint
-JOB2LOAD= 01;
+JOB2LOAD= 0;
 %% OPTIONS
 SIMID   = 'ZP';  % Name of the simulation
 CO      = -2;  % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)