diff --git a/matlab/COSOlver_matrices_analysis.m b/matlab/COSOlver_matrices_analysis.m
new file mode 100644
index 0000000..1514d48
--- /dev/null
+++ b/matlab/COSOlver_matrices_analysis.m
@@ -0,0 +1,72 @@
+addpath(genpath('../matlab')) % ... add
+%% Grid configuration
+N       = 10;     % Frequency gridpoints (Nkr = N/2)
+L       = 120;     % Size of the squared frequency domain
+dk      = 2*pi/L;
+kmax    = N/2*dk;
+kr      = dk*(0:N/2);
+kz      = dk*(0:N/2);
+[KZ, KR]= meshgrid(kz,kr);
+KPERP   = sqrt(KR.^2 + KZ.^2);
+kperp   = reshape(KPERP,[1,numel(kr)^2]);
+kperp   = uniquetol(kperp,1e-14);
+Nperp   = numel(kperp);
+COSOlver_path = '../../Documents/MoliSolver/COSOlver/';
+COSOLVER.pmaxe = 10;
+COSOLVER.jmaxe = 5;
+COSOLVER.pmaxi = 10;
+COSOLVER.jmaxi = 5;
+
+COSOLVER.neFLR  = max(5,ceil(COSOLVER.kperp^2)); % rule of thumb for sum truncation
+COSOLVER.niFLR  = max(5,ceil(COSOLVER.kperp^2));
+
+COSOLVER.neFLRs = 0; %  ... only for GK abel 
+COSOLVER.npeFLR = 0; %  ... only for GK abel 
+COSOLVER.niFLRs = 0; %  ... only for GK abel 
+COSOLVER.npiFLR = 0; %  ... only for GK abel 
+
+COSOLVER.gk = 1; 
+COSOLVER.co = 3;
+if 0
+    %% plot the kperp distribution
+   figure
+   plot(kperp)
+end
+%% Load the matrices
+C_self_i = zeros((COSOLVER.pmaxi+1)*(COSOLVER.jmaxi+1),(COSOLVER.pmaxi+1)*(COSOLVER.jmaxi+1),Nperp);
+
+for n_ = 1:Nperp
+    COSOLVER.kperp = kperp(n_);
+    k_string      = sprintf('%0.4f',kperp(n_));
+    mat_file_name = ['../iCa/self_Coll_GKE_',num2str(COSOLVER.gk),'_GKI_',num2str(COSOLVER.gk),...
+    '_ESELF_',num2str(COSOLVER.co),'_ISELF_',num2str(COSOLVER.co),...
+    '_Pmaxe_',num2str(COSOLVER.pmaxe),'_Jmaxe_',num2str(COSOLVER.jmaxe),...
+    '_Pmaxi_',num2str(COSOLVER.pmaxi),'_Jmaxi_',num2str(COSOLVER.jmaxi),...
+    '_JE_12',...
+    '_NFLR_',num2str(COSOLVER.neFLR),...
+    '_kperp_',k_string,'.h5'];
+
+    tmp  = h5read(mat_file_name,'/Caapj/Ceepj');
+    C_self_i(:,:,n_) = tmp;
+end
+
+%% Post processing
+dC_self_i = diff(C_self_i,1,3);
+gvar_dC    = squeeze(sum(sum(abs(dC_self_i),2),1));
+%% Plots
+%% all coeffs evolution
+figure
+for ip_ = 1:COSOLVER.pmaxi+1
+    for ij_ = 1:COSOLVER.jmaxi+1
+        plot(kperp,squeeze(C_self_i(ip_,ij_,:)),'o'); hold on;
+    end
+end
+%% global matrix variation
+figure;
+kperpmid = 0.5*(kperp(2:end)+kperp(1:end-1));
+plot(kperpmid,gvar_dC./diff(kperp)');
+
+
+
+
+
diff --git a/matlab/write_fort90_COSOlver.m b/matlab/write_fort90_COSOlver.m
new file mode 100644
index 0000000..46da401
--- /dev/null
+++ b/matlab/write_fort90_COSOlver.m
@@ -0,0 +1,84 @@
+INPUT = [COSOlver_path,'fort.90'];
+fid = fopen(INPUT,'wt');
+
+fprintf(fid,'&BASIC\n');
+fprintf(fid,['pmaxe = ', num2str(COSOLVER.pmaxe),'\n']);
+fprintf(fid,['jmaxe = ', num2str(COSOLVER.jmaxe),'\n']);
+fprintf(fid,['pmaxi = ', num2str(COSOLVER.pmaxi),'\n']);
+fprintf(fid,['jmaxi = ', num2str(COSOLVER.jmaxi),'\n']);
+fprintf(fid,'JEmaxx=12\n');
+fprintf(fid,'PMmaxx=16\n');
+fprintf(fid,'\n');
+
+fprintf(fid,['neFLR=',num2str(COSOLVER.neFLR),'\n']);
+fprintf(fid,['neFLRs =',num2str(COSOLVER.neFLRs),'\n']);
+fprintf(fid,['npeFLR=',num2str(COSOLVER.npeFLR),'\n']);
+fprintf(fid,['niFLR=',num2str(COSOLVER.niFLR),'\n']);
+fprintf(fid,['niFLRs=',num2str(COSOLVER.niFLRs),'\n']);
+fprintf(fid,['npiFLR=',num2str(COSOLVER.npiFLR),'\n']);
+fprintf(fid,'\n');
+fprintf(fid,['eecolls=.true.','\n']);
+fprintf(fid,['iicolls=.true.','\n']);
+fprintf(fid,['eicolls=.true.','\n']);
+fprintf(fid,['iecolls=.true.','\n']);
+fprintf(fid,'/\n');
+fprintf(fid,'\n');
+
+fprintf(fid,'&BASIS_TRANSFORMATION_PAR\n');
+fprintf(fid,['T5dir = ','''/misc/coeffs_backup/T5src/''','\n']);
+fprintf(fid,['T4dir = ','''/misc/T4/NNT4_L000x200_K000x200_P000x200_J000x200/''','\n']);
+fprintf(fid,'idxT4max = 30\n');
+fprintf(fid,'idxT5max = 0\n');
+fprintf(fid,'IFT4 = .true.\n');
+fprintf(fid,'IFT5 = .false.\n');
+fprintf(fid,'/\n');
+fprintf(fid,'\n');
+
+fprintf(fid,'&MODEL_PAR \n');
+fprintf(fid,'nu=1\n');
+fprintf(fid,'mu=0.023338\n');
+fprintf(fid,'tau=1\n');
+fprintf(fid,['kperp=',num2str(COSOLVER.kperp,16),'\n']);
+fprintf(fid,'/\n');
+fprintf(fid,'\n');
+
+fprintf(fid,'&OPERATOR_MODEL \n');
+fprintf(fid,['ETEST=', num2str(COSOLVER.co),'\n']);
+fprintf(fid,['EBACK=', num2str(COSOLVER.co),'\n']);
+fprintf(fid,['ITEST=', num2str(COSOLVER.co),'\n']);
+fprintf(fid,['IBACK=', num2str(COSOLVER.co),'\n']);
+fprintf(fid,'\n');
+fprintf(fid,['ESELF=', num2str(COSOLVER.co),'\n']);
+fprintf(fid,['ISELF=', num2str(COSOLVER.co),'\n']);
+fprintf(fid,'\n');
+fprintf(fid,['GKE = ', num2str(COSOLVER.gk),'\n']);
+fprintf(fid,['GKI = ', num2str(COSOLVER.gk),'\n']);
+fprintf(fid,'DKTEST = .F.\n');
+fprintf(fid,'DKBACK = .F.\n');
+fprintf(fid,'ADDTEST = .T.\n');
+fprintf(fid,'ADDBACK = .T.\n');
+fprintf(fid,'only_gk_part = .false.\n');
+fprintf(fid,'only_symmetric = .false.\n');
+fprintf(fid,'/\n');
+fprintf(fid,'\n');
+
+fprintf(fid,'&MPI\n');
+fprintf(fid,'Nprocs_j_ii = 1\n');
+fprintf(fid,'MPI_balanced = .false.\n');
+fprintf(fid,'/\n');
+fprintf(fid,'\n');
+
+fprintf(fid,'&OUTPUT_PAR\n');
+fprintf(fid,'OVERWRITE=.true.\n');
+fprintf(fid,'nsave_ei=0\n');
+fprintf(fid,'nsave_ie=0\n');
+fprintf(fid,'nsave_ee=0\n');
+fprintf(fid,'nsave_ii=2\n');
+fprintf(fid,'ifrestart=.false.\n');
+fprintf(fid,['suffix_resfile=','''''','\n']);
+fprintf(fid,['outdir = ','''../../../HeLaZ/iCa''','\n']);
+fprintf(fid,'/\n');
+fprintf(fid,'\n');
+
+
+fclose(fid);
\ No newline at end of file
diff --git a/wk/compute_collision_mat.m b/wk/compute_collision_mat.m
new file mode 100644
index 0000000..4696074
--- /dev/null
+++ b/wk/compute_collision_mat.m
@@ -0,0 +1,81 @@
+addpath(genpath('../matlab')) % ... add
+%% Grid configuration
+N       = 200;     % Frequency gridpoints (Nkr = N/2)
+L       = 120;     % Size of the squared frequency domain
+dk      = 2*pi/L;
+kmax    = N/2*dk;
+kr      = dk*(0:N/2);
+kz      = dk*(0:N/2);
+[KZ, KR]= meshgrid(kz,kr);
+KPERP   = sqrt(KR.^2 + KZ.^2);
+kperp   = reshape(KPERP,[1,numel(kr)^2]);
+kperp   = uniquetol(kperp,1e-14);
+Nperp   = numel(kperp);
+if 0
+    %% plot the kperp distribution
+   figure
+   plot(kperp)
+end
+%% Check if the differences btw kperp is larger than naming precision
+dkperp  = diff(kperp);
+warning = sum(dkperp<0.0002);
+if warning > 0
+    disp('Warning : dkperp < 0.0002');
+end
+%%
+n_ = 1;
+for k_ = kperp
+    %% Script to run COSOlver in order to create needed collision matrices
+    COSOlver_path = '../../Documents/MoliSolver/COSOlver/';
+    COSOLVER.pmaxe = 10;
+    COSOLVER.jmaxe = 5;
+    COSOLVER.pmaxi = 10;
+    COSOLVER.jmaxi = 5;
+    COSOLVER.kperp = k_;
+
+    COSOLVER.neFLR  = max(5,ceil(COSOLVER.kperp^2)); % rule of thumb for sum truncation
+    COSOLVER.niFLR  = max(5,ceil(COSOLVER.kperp^2));
+
+    COSOLVER.neFLRs = 0; %  ... only for GK abel 
+    COSOLVER.npeFLR = 0; %  ... only for GK abel 
+    COSOLVER.niFLRs = 0; %  ... only for GK abel 
+    COSOLVER.npiFLR = 0; %  ... only for GK abel 
+
+    COSOLVER.gk = 1; 
+    COSOLVER.co = 3;
+    % ! 0 = nothing 
+    % ! 1 = coulomb
+    % ! 2 = pitch-angle (with constant coll.freq.)
+    % ! 3 = sugama
+    % ! 4 = pitch-angle with veloty dependent coll. freq.
+    % ! 5 = improved sugama
+    % ! 6 = Hirschman-Sigmar Clarke
+    % ! 7 = Abel (only for like species)
+    % ! 8 = conserving momentun pitch-angle operator (with veloty dependent coll. freq.)
+    % ! 9 = GK coulomb polarization term
+
+    write_fort90_COSOlver
+    
+    k_string      = sprintf('%0.4f',k_);
+    mat_file_name = ['../iCa/self_Coll_GKE_',num2str(COSOLVER.gk),'_GKI_',num2str(COSOLVER.gk),...
+        '_ESELF_',num2str(COSOLVER.co),'_ISELF_',num2str(COSOLVER.co),...
+        '_Pmaxe_',num2str(COSOLVER.pmaxe),'_Jmaxe_',num2str(COSOLVER.jmaxe),...
+        '_Pmaxi_',num2str(COSOLVER.pmaxi),'_Jmaxi_',num2str(COSOLVER.jmaxi),...
+        '_JE_12',...
+        '_NFLR_',num2str(COSOLVER.neFLR),...
+        '_kperp_',k_string,'.h5'];
+    if (exist(mat_file_name,'file') > 0)
+        disp(['Matrix available for kperp = ',k_string]);
+    else
+        cd ../../Documents/MoliSolver/COSOlver/
+        disp(['Matrix not found for kperp = ',k_string]);
+        disp([num2str(n_),'/',Nperp]
+        disp('computing...');
+        CMD = 'mpirun -np 6 bin/CO 2 2 2 > out.txt';
+        disp(CMD); 
+        system(CMD);
+        system(CMD);
+        disp('..done');
+        cd ../../../HeLaZ/wk
+    end
+end
\ No newline at end of file