diff --git a/matlab/default_plots_options.m b/matlab/default_plots_options.m
new file mode 100644
index 0000000000000000000000000000000000000000..c4cd384709401599ada43bcad1f7570cf9241380
--- /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 0000000000000000000000000000000000000000..550679e27c754b6770c0f4f2181d9eeaf6d57642
--- /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 0000000000000000000000000000000000000000..308f58d6b8f73ad4ced95301d059704904e9ec62
--- /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