diff --git a/matlab/default_plots_options.m b/matlab/default_plots_options.m
new file mode 100644
index 0000000..c4cd384
--- /dev/null
+++ b/matlab/default_plots_options.m
@@ -0,0 +1,9 @@
+%% Default values for plots
+set(0,'defaulttextInterpreter','latex') 
+set(0, 'defaultAxesTickLabelInterpreter','latex');
+set(0, 'defaultLegendInterpreter','latex');
+set(0,'defaultAxesFontSize',16)
+set(0, 'DefaultLineLineWidth', 1.5);
+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
diff --git a/matlab/helaz_analysis.m b/matlab/helaz_analysis.m
new file mode 100644
index 0000000..550679e
--- /dev/null
+++ b/matlab/helaz_analysis.m
@@ -0,0 +1,112 @@
+%% HeLaZ data
+filename = 'results_00.h5';
+default_plots_options % Script to set up default plot variables
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Load the data
+moment = 'Ni00';
+
+kr       = h5read(filename,['/data/var2d/' moment '/coordkr']);
+kz       = h5read(filename,['/data/var2d/' moment '/coordkz']);
+timeNi   = h5read(filename,'/data/var2d/time');
+Nipj     = zeros(numel(timeNi),numel(kr),numel(kz));
+for it = 1:numel(timeNi)
+    tmp          = h5read(filename,['/data/var2d/', moment,'/', num2str(it,'%06d')]);
+    Nipj(it,:,:) = tmp.real + 1i * tmp.imaginary; 
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Plot growth rate vs t
+gammas = zeros(numel(kr),numel(kz));
+shifts = zeros(numel(kr),numel(kz));
+% Linear fit of log(Napj)
+x1    = timeNi;
+itmin = ceil(0.5 * numel(timeNi)); %Take the second half of the time evolution
+for ikr = 1:numel(kr)
+    for ikz = 1:numel(kz)
+        fit = polyfit(x1(itmin:end),log(abs(Nipj(itmin:end,ikr,ikz))),1);
+        gammas(ikr,ikz) = fit(1);
+        shifts(ikr,ikz) = fit(2);
+    end
+end
+
+FIGNAME = 'gamma_t';
+fig = figure;
+for ikr = 1:numel(kr)
+    linename = ['$k_r = ',num2str(kr(ikr)),'$'];
+    plot(kz,gammas(ikr,:),'DisplayName',linename);
+end
+TITLE  = [];
+TITLE = [TITLE,'$\eta_n=',num2str(1.0/MODEL.eta_n),'$, '];
+TITLE = [TITLE,'$\eta_B=',num2str(MODEL.eta_B),'$, '];
+TITLE = [TITLE,   '$\nu=',num2str(MODEL.nu),'$, '];
+%TITLE = [TITLE,   '$k_z=',num2str(GRID.kz),'$'];
+
+title(TITLE);
+grid on
+legend('show')
+xlabel('$t$')
+ylabel(['$|',moment,'|$'])
+
+%% Saving fig
+if SAVEFIG
+    FIGDIR = ['../results/', SIMID,'/'];
+    if ~exist(FIGDIR, 'dir')
+       mkdir(FIGDIR)
+    end
+    FIGNAME = [FIGNAME,'_Pe_',num2str(GRID.pmaxe),'_Je_',num2str(GRID.jmaxe),...
+        '_Pi_',num2str(GRID.pmaxi),'_Ji_',num2str(GRID.jmaxi),...
+        '_etan_',num2str(MODEL.eta_n),'_etaB_',num2str(MODEL.eta_B),'_nu_',num2str(MODEL.nu)];
+    FIGNAME = [FIGDIR, FIGNAME,'.fig'];
+    savefig(fig,FIGNAME);
+    disp(['Figure saved @ : ',FIGNAME])
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Plot Ni00 evolution
+FIGNAME = 'Ni00_t';
+fig = figure;
+%HeLaZ results
+x1 = timeNi;
+il = 1;
+for ikr = 1:numel(kr)
+    ic = 1;
+    for ikz = 1:2:numel(kz)
+        linename = ['$k_r = ',num2str(kr(ikr)),', k_z = ', num2str(kz(ikz)),'$'];
+        y1 = abs(Nipj(:,ikr,ikz));
+        semilogy(x1,y1,...
+            'DisplayName',linename,'Color', line_colors(ic,:), 'LineStyle', '--')
+         hold on
+        semilogy(x1(itmin:end),exp(gammas(ikr,ikz)*x1(itmin:end) + shifts(ikr,ikz)),...
+            'DisplayName',[linename,' fit'],'Color', line_colors(ic,:), 'LineStyle', '-.')
+        ic = ic + 1;
+    end
+    il = il + 1;
+end
+
+TITLE  = [];
+TITLE = [TITLE,'$\eta_n=',num2str(1.0/MODEL.eta_n),'$, '];
+TITLE = [TITLE,'$\eta_B=',num2str(MODEL.eta_B),'$, '];
+TITLE = [TITLE,   '$\nu=',num2str(MODEL.nu),'$, '];
+%TITLE = [TITLE,   '$k_z=',num2str(GRID.kz),'$'];
+
+title(TITLE);
+grid on
+legend('show')
+xlabel('$t$')
+ylabel(['$|',moment,'|$'])
+
+%% Saving fig
+if SAVEFIG
+    FIGDIR = ['../results/', SIMID,'/'];
+    if ~exist(FIGDIR, 'dir')
+       mkdir(FIGDIR)
+    end
+    FIGNAME = [FIGNAME,'_Pe_',num2str(GRID.pmaxe),'_Je_',num2str(GRID.jmaxe),...
+        '_Pi_',num2str(GRID.pmaxi),'_Ji_',num2str(GRID.jmaxi),...
+        '_etan_',num2str(MODEL.eta_n),'_etaB_',num2str(MODEL.eta_B),'_nu_',num2str(MODEL.nu)];
+    FIGNAME = [FIGDIR, FIGNAME,'.fig'];
+    savefig(fig,FIGNAME);
+    disp(['Figure saved @ : ',FIGNAME])
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\ No newline at end of file
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
new file mode 100644
index 0000000..308f58d
--- /dev/null
+++ b/matlab/write_fort90.m
@@ -0,0 +1,58 @@
+function [INPUT] = write_fort90(OUTPUTS,GRID,MODEL,INITIAL,TIME_INTEGRATION,BASIC)
+% Write the input script "fort.90" with desired parameters
+INPUT = '../wk/fort.90';
+fid = fopen(INPUT,'wt');
+fprintf(fid,'&BASIC\n');
+fprintf(fid,['  nrun=', num2str(BASIC.nrun),'\n']);
+fprintf(fid,['  dt=', num2str(BASIC.dt),'\n']);
+fprintf(fid,['  tmax=', num2str(BASIC.tmax),'    ! time normalized to 1/omega_pe\n']);
+fprintf(fid,'/\n');
+fprintf(fid,'&GRID\n');
+fprintf(fid,['  pmaxe =', num2str(GRID.pmaxe),'    ! number of Hermite moments \n']);
+fprintf(fid,['  jmaxe = ', num2str(GRID.jmaxe),'     ! number of Hermite moments \n']);
+fprintf(fid,['  pmaxi = ', num2str(GRID.pmaxi),'    ! number of Hermite moments \n']);
+fprintf(fid,['  jmaxi = ', num2str(GRID.jmaxi),'     ! number of Hermite moments\n']);
+fprintf(fid,['  nkr   = ', num2str(GRID.nkr),'\n']);
+fprintf(fid,['  krmin = ', num2str(GRID.krmin),'\n']);
+fprintf(fid,['  krmax = ', num2str(GRID.krmax),'     ! Normalized to cs0/Omega_i\n']);
+fprintf(fid,['  nkz   = ', num2str(GRID.nkz),'\n']);
+fprintf(fid,['  kzmin = ', num2str(GRID.kzmin),'\n']);
+fprintf(fid,['  kzmax = ', num2str(GRID.kzmax),'     ! Normalized to cs0/Omega_i\n']);
+fprintf(fid,'/\n');
+fprintf(fid,'&OUTPUT_PAR\n');
+fprintf(fid,['  nsave_0d = ', num2str(OUTPUTS.nsave_0d),'\n']);
+fprintf(fid,['  nsave_1d = ', num2str(OUTPUTS.nsave_1d),'\n']);
+fprintf(fid,['  nsave_2d = ', num2str(OUTPUTS.nsave_2d),'\n']);
+fprintf(fid,['  nsave_5d = ', num2str(OUTPUTS.nsave_5d),'\n']);
+fprintf(fid,['  write_Ni00    = ', OUTPUTS.write_Ni00,'\n']);
+fprintf(fid,['  write_moments = ', OUTPUTS.write_moments,'\n']);
+fprintf(fid,['  write_phi     = ', OUTPUTS.write_phi,'\n']);
+fprintf(fid,['  write_doubleprecision = ', OUTPUTS.write_doubleprecision,'\n']);
+fprintf(fid,['  resfile0      = ', OUTPUTS.resfile0,'\n']);
+fprintf(fid,'/\n');
+fprintf(fid,'&MODEL_PAR\n');
+fprintf(fid,'  ! Collisionality\n');
+fprintf(fid,['  nu      = ', num2str(MODEL.nu),'       ! Normalized collision frequency normalized to plasma frequency\n']);
+fprintf(fid,['  tau_e   = ', num2str(MODEL.tau_e),'         ! T_e/T_e\n']);
+fprintf(fid,['  tau_i   = ', num2str(MODEL.tau_i),'         ! T_i/T_e       temperature ratio\n']);
+fprintf(fid,['  sigma_e = ', num2str(MODEL.sigma_e),'   ! sqrt(m_e/m_i) mass ratio\n']);
+fprintf(fid,['  sigma_i = ', num2str(MODEL.sigma_i),'         ! sqrt(m_i/m_i)\n']);
+fprintf(fid,['  q_e     = ', num2str(MODEL.q_e),'         ! Electrons charge\n']);
+fprintf(fid,['  q_i     = ', num2str(MODEL.q_i),'         ! Ions charge\n']);
+fprintf(fid,['  eta_n   = ', num2str(MODEL.eta_n),'         ! L_perp / L_n Density gradient\n']);
+fprintf(fid,['  eta_T   = ', num2str(MODEL.eta_T),'         ! L_perp / L_T Temperature gradient\n']);
+fprintf(fid,['  eta_B   = ', num2str(MODEL.eta_B),'         ! L_perp / L_B Magnetic gradient and curvature\n']);
+fprintf(fid,['  lambdaD = ', num2str(MODEL.lambdaD),'         ! Debye length\n']);
+fprintf(fid,'/\n');
+fprintf(fid,'&INITIAL_CON\n');
+fprintf(fid,'  ! Background value\n');
+fprintf(fid,['  initback_moments=', num2str(INITIAL.initback_moments),' ! x 1e-3\n']);
+fprintf(fid,'  ! Noise amplitude\n');
+fprintf(fid,['  initnoise_moments=', num2str(INITIAL.initnoise_moments),'\n']);
+fprintf(fid,['  iseed=', num2str(INITIAL.iseed),'\n']);
+fprintf(fid,'/\n');
+fprintf(fid,'&TIME_INTEGRATION_PAR\n');
+fprintf(fid,['  numerical_scheme=', TIME_INTEGRATION.numerical_scheme,'\n']);
+fprintf(fid,'/');
+fclose(fid);
+end