diff --git a/quad10/fenomens/lab/p4/MC-2.f90 b/quad10/fenomens/lab/p4/MC-2.f90
index 554890d..ccc2941 100644
--- a/quad10/fenomens/lab/p4/MC-2.f90
+++ b/quad10/fenomens/lab/p4/MC-2.f90
@@ -7,19 +7,25 @@
   integer*4 :: I, J, IMC, MCTOT, MCINI, MCD, IPAS, N, SEED, NSEED, SEED0, pbc
   integer :: genrand_int31
   character(100) :: NOM
-  logical :: SAVEFINALCONF
+  logical :: SAVEFINALCONF, SAVEEVOLUTION
 
   integer :: SUMI
-  real(dp) :: SUME, SUME2, SUMM, SUMAM, SUMM2, VARE, VARM, CV, CHI
+  real(dp) :: SUME, SUME2, SUMM, SUMAM, SUMM2, VARE_TIMES_SUMI, VARM_TIMES_SUMI, CV, CHI
 
-  ! This should be replaced by L but it doesn't work dynamically, so we're
-  ! setting 64 as the upper limit which should be higher than L.
-  integer*2 :: S(1:64, 1:64)
+  ! Donat que les matrius a Fortran han de tenir una mida fixada i no poden
+  ! dependre dinàmicament d'un paràmetre L, fem S de mida 128x128, i establim
+  ! 128 com el valor màxim de L.
+  integer*2 :: S(1:128, 1:128)
 
-  ! Inicialitzem algunes variables sobre el problema
-  namelist /DADES/ L, NOM, TEMP, NSEED, SEED0, MCTOT, MCINI, MCD, SAVEFINALCONF
+  ! Inicialitzem algunes variables sobre el problema, que es passen al programa amb un namelist
+  namelist /DADES/ L, NOM, TEMP, NSEED, SEED0, MCTOT, MCINI, MCD, SAVEFINALCONF, SAVEEVOLUTION
   read(nml = DADES, unit = 5)
 
+  if (L < 1 .or. L > 128) then
+    print *, "ERROR: L ha de ser un enter entre 1 i 128."
+    stop
+  endif
+
   N = L*L
 
   ! Cache dels valors de l'exponencial
@@ -27,6 +33,7 @@
     W(I) = (2**30 - 1 + 2**30)*exp(-float(I)/TEMP)
   enddo
 
+  ! Variables on anirem sumant les magnituds per fer el promig posteriorment
   SUMI = 0
   SUME = 0
   SUME2 = 0
@@ -35,6 +42,12 @@
   SUMAM = 0
   SUMM2 = 0
 
+  ! Si s'ha especificat de guardar l'evolució temporal, obrim el fitxer on la guardarem
+  if (SAVEEVOLUTION) then
+    open(12, file = "data_out/" // trim(NOM) // ".ev")
+  endif
+
+  ! Fem la simulació per NSEED seeds inicials.
   do SEED = SEED0, SEED0 + NSEED - 1
     print *, "Seed: ", SEED
 
@@ -43,11 +56,13 @@
 
     E = energ(S, L)
     M = magne(S, L)
+    if (SAVEEVOLUTION) then
+      write (12, *) "0", E, M
+    endif
 
     ! Iterem amb el mètode de Montecarlo
     do IMC = 1, MCTOT
       do IPAS = 1, N
-        ! LOS LOOPS HACEN COSAS
         I = mod(genrand_int31(), L) + 1
         J = mod(genrand_int31(), L) + 1
         SUMA = S(PBC(I - 1, L), J) + &
@@ -59,7 +74,7 @@
         if (DIFE > 0) then
           DELTA = genrand_int31()
           if (DELTA >= W(int(DIFE))) then
-            ! NO ho acceptem
+            ! NO ho acceptem, passem a la següent volta del loop
             cycle
           endif
         endif
@@ -70,6 +85,12 @@
         M = M + 2*S(I, J)
       enddo
 
+      if (SAVEEVOLUTION) then
+        write (12, *) IMC, E, M
+      endif
+
+      ! Per evitar les correlacions de configuracions consecutives, només
+      ! prenem mesures cada MCD configuracions.
       if (IMC > MCINI .and. mod(IMC, MCD) == 0) then
         SUMI = SUMI + 1
 
@@ -89,15 +110,15 @@
   SUMM = SUMM/SUMI
   SUMAM = SUMAM/SUMI
   SUMM2 = SUMM2/SUMI
-  VARE = SUME2 - SUME**2
-  VARM = SUMM2 - SUMM**2
+  VARE_TIMES_SUMI = SUME2 - SUME**2
+  VARM_TIMES_SUMI = SUMM2 - SUMM**2
   CV = (SUME2 - SUME**2)/real(N*TEMP**2)
   CHI = (SUMM2 - SUMAM**2)/real(N*TEMP)
 
-  ! Guardem a un fiter els promitjos
+  ! Guardem a un fitxer els promitjos
   open(unit = 13, file = "data_out/" // trim(NOM) // ".res")
-  ! Valors: L, T, <E>, <E**2>, Var(E), <M>, <|M|>, <M**2>, Var(M), C_V, CHI
-  write(13, *) L, TEMP, SUME, SUME2, VARE, SUMM, SUMAM, SUMM2, VARM, CV, CHI
+  ! Valors: L, T, <E>, <E**2>, Var(E)*SUMI, <M>, <|M|>, <M**2>, Var(M)*SUMI, C_V, CHI
+  write(13, *) L, TEMP, SUME, SUME2, VARE_TIMES_SUMI, SUMM, SUMAM, SUMM2, VARM_TIMES_SUMI, CV, CHI
 
   if (SAVEFINALCONF) then
     call writeConfig(S, L, "data_out/" // trim(NOM) // ".conf")
@@ -105,4 +126,7 @@
 
   ! Tanquem els fitxers de sortida
   close(13)
+  if (SAVEEVOLUTION) then
+    close(12)
+  endif
 endprogram main
diff --git a/quad10/fenomens/lab/p4/Makefile b/quad10/fenomens/lab/p4/Makefile
index 21e5cbf..9c7b58a 100644
--- a/quad10/fenomens/lab/p4/Makefile
+++ b/quad10/fenomens/lab/p4/Makefile
@@ -2,9 +2,11 @@
 FC=gfortran
 FC_FLAGS=-O3 -Wall #-g -fbacktrace -fcheck=all
 
-all: MC-2.o
+all: MC-2.o out/cv_via_derivative
 make_out_folder:
 	mkdir -p out data_out
+out/cv_via_derivative: cv_via_derivative.cc
+	g++ cv_via_derivative.cc -Wall -o out/cv_via_derivative
 MC-2.o: make_out_folder mt19937ar.o generateSpinMatrix.o magne.o energ.o pbc.o writeconfig.o MC-2.f90
 	$(FC) $(FC_FLAGS) -c MC-2.f90
 	$(FC) $(FC_FLAGS) MC-2.o generateSpinMatrix.o mt19937ar.o magne.o energ.o pbc.o writeconfig.o -o out/mc2
diff --git a/quad10/fenomens/lab/p4/README.md b/quad10/fenomens/lab/p4/README.md
index 421cd3f..47b0ec3 100644
--- a/quad10/fenomens/lab/p4/README.md
+++ b/quad10/fenomens/lab/p4/README.md
@@ -2,8 +2,21 @@
 **Data: 12 de maig de 2022**
 
 ## Instruccions per compilar
-Per compilar, executeu `make all`. Els executables es trobaran a la carpeta `out`.
+Per compilar, executeu `make all`. L'executable de la simulació es trobarà a la carpeta `out`, tot i que no s'hauria de cridar aquest executable directament.
 
 ## Programes
-- `run.bash`: executa la simulació per un valor de T i L concrets. Per més informació, executeu `./run.bash --help`.
+*** note
+**Atenció!** Per poder executar aquests programes, s'ha d'haver compilat abans l'executable del codi Fortran com s'especifica a la secció anterior.
+***
+
+Per tal d'executar la simulació, existeixen diversos scripts de bash que es poden cridar:
+- `run.bash`: executa la simulació per un valor de T i L concrets, amb un munt de valors addicionals configurables (nombre de llavors, nombre d'iteracions, etc.). Per més informació, executeu `./run.bash --help`.
 - `runAll.bash`: executa les simulacions per diferents valors de T i L, i guarda les dades finals a la carpeta `data_out/dep_en_L`.
+
+També existeixen diversos scipts per tal de generar algunes gràfiques (altes gràfiques es generen directament amb scripts que hi ha a la [carpeta de l'informe](../informe/)):
+- `drawDepEnLGraphs.bash`: dibuixa les gràfiques de la dependència de certes variables amb la temperatura per diferents Ls.
+   - Precondició: haver executat `runAll.bash` abans, tot i que sigui parcialment.
+- `drawFinalConfsGraph.bash`: dibuixa un diagrama amb les configuracions finals per diferents valors de T.
+   - No té cap precondició.
+- `drawDepEnSeedGraph.bash`: dibuixa un gràfic que mostra la dependència de l'energia i la magnetització amb l'elecció de la seed inicial.
+   - No té cap precondició.
diff --git a/quad10/fenomens/lab/p4/cv_via_derivative.cc b/quad10/fenomens/lab/p4/cv_via_derivative.cc
new file mode 100644
index 0000000..c5308a2
--- /dev/null
+++ b/quad10/fenomens/lab/p4/cv_via_derivative.cc
@@ -0,0 +1,29 @@
+#include <iostream>
+#include <vector>
+using namespace std;
+
+// Fitxer que calcula a partir d'un fitxer de resultats del programa mc-2 els
+// valors de d<e>/dT amb derivació numèrica.
+
+struct Row {
+  long double T;
+  long double E;
+};
+
+int main() {
+  vector<Row> data;
+  long double t; // Trash
+  long double T;
+  long double E;
+  while (cin >> t >> T >> E >> t >> t >> t >> t >> t >> t >> t >> t) {
+    Row r;
+    r.T = T;
+    r.E = E;
+    data.push_back(r);
+  }
+
+  for (int i = 0; i < (int)data.size() - 1; ++i) {
+    long double derivada = (data[i + 1].E - data[i].E)/(data[i + 1].T - data[i].T);
+    cout << data[i].T << " " << derivada << endl;
+  }
+}
diff --git a/quad10/fenomens/lab/p4/drawDepEnLGraphs.bash b/quad10/fenomens/lab/p4/drawDepEnLGraphs.bash
new file mode 100755
index 0000000..08914a8
--- /dev/null
+++ b/quad10/fenomens/lab/p4/drawDepEnLGraphs.bash
@@ -0,0 +1,28 @@
+#!/bin/bash
+cd data_out/dep_en_L
+LList="8 16 32 64 96"
+LFinalList=()
+
+mkdir -p tmpdata
+for L in $LList; do
+  if [ -d "$L" ]; then
+    (cd $L; cat *.res | sort -nk2 > "../tmpdata/$L.dat")
+    if [ $L -eq 32 ]; then
+      ../../out/cv_via_derivative < tmpdata/32.dat > tmpdata/32_derivada_E.dat
+    fi
+    LFinalList+=("$L")
+  else
+    echo "Folder with L=$L doesn't exist, skipping."
+  fi
+done
+
+LFinalListString="${LFinalList[@]}"
+
+cd ../../
+mkdir -p data_out/dep_en_L_graphs
+
+cd data_out/dep_en_L_graphs
+../../graphs/dependenciaEnL.gnu "$LFinalListString"
+
+cd ../../
+#rm -rf data_out/dep_en_L/tmpdata
diff --git a/quad10/fenomens/lab/p4/drawDepEnSeedGraph.bash b/quad10/fenomens/lab/p4/drawDepEnSeedGraph.bash
new file mode 100755
index 0000000..cf1a6ab
--- /dev/null
+++ b/quad10/fenomens/lab/p4/drawDepEnSeedGraph.bash
@@ -0,0 +1,7 @@
+#!/bin/bash
+T=1.9
+L=64
+SEEDS="17427 9090712 590100 8248 99963"
+mkdir -p data_out/depEnSeed
+parallel -P -1 --linebuffer --tag "./run.bash $T $L depEnSeed/{} --saveEvolution --mcTot 10000 --nSeeds 1 --initialSeed {} --skipIfComputed" ::: $SEEDS
+(cd data_out/depEnSeed; ../../graphs/dependenciaEnSeed.gnu "$SEEDS" "$L")
diff --git a/quad10/fenomens/lab/p4/drawFinalConfsGraph.bash b/quad10/fenomens/lab/p4/drawFinalConfsGraph.bash
new file mode 100755
index 0000000..7f83e9d
--- /dev/null
+++ b/quad10/fenomens/lab/p4/drawFinalConfsGraph.bash
@@ -0,0 +1,6 @@
+#!/bin/bash
+TList=$(seq 1.9 0.3 3.4)
+L=64
+mkdir -p data_out/finalconfs
+parallel -P -1 --linebuffer --tag "./run.bash {} $L finalconfs/{} --saveFinalConf --nSeeds 1 --initialSeed 99963 --skipIfComputed" ::: $TList
+(cd data_out/finalconfs; ../../graphs/finalConfs.gnu "$TList" "$L")
diff --git a/quad10/fenomens/lab/p4/drawGraphs.bash b/quad10/fenomens/lab/p4/drawGraphs.bash
deleted file mode 100755
index 9a55c5c..0000000
--- a/quad10/fenomens/lab/p4/drawGraphs.bash
+++ /dev/null
@@ -1,21 +0,0 @@
-#!/bin/bash
-cd data_out/dep_en_L
-LList="8 16 32 64"
-LFinalList=()
-
-mkdir -p tmpdata
-for L in $LList; do
-  if [ -d "$L" ]; then
-    (cd $L; cat *.res | sort -nk2 > "../tmpdata/$L.dat")
-    LFinalList+=("$L")
-  else
-    echo "Folder with L=$L doesn't exist, skipping."
-  fi
-done
-
-LFinalListString="${LFinalList[@]}"
-
-cd ../../
-mkdir -p data_out/dep_en_L_graphs
-./graphs/dependenciaEnL.gnu "$LFinalListString"
-rm -rf data_out/dep_en_L/tmpdata
diff --git a/quad10/fenomens/lab/p4/energ.f90 b/quad10/fenomens/lab/p4/energ.f90
index 431a296..a57601e 100644
--- a/quad10/fenomens/lab/p4/energ.f90
+++ b/quad10/fenomens/lab/p4/energ.f90
@@ -1,5 +1,7 @@
+! Retorna, donada una matriu d'spins S d'un sistema de mida LxL,
+! l'energia del sistema.
 real*8 function energ(S, L)
-  integer*2 :: S(1:64, 1:64)
+  integer*2 :: S(1:128, 1:128)
   integer*4 :: I, J, L, PBC
   real*8 :: ENE
   ENE = 0.0d0
diff --git a/quad10/fenomens/lab/p4/generateSpinMatrix.f90 b/quad10/fenomens/lab/p4/generateSpinMatrix.f90
index 0128119..6351e45 100644
--- a/quad10/fenomens/lab/p4/generateSpinMatrix.f90
+++ b/quad10/fenomens/lab/p4/generateSpinMatrix.f90
@@ -1,9 +1,12 @@
+! Inicialitza la matriu d'spins S aleatòriament, donada la mida del sistema L
+! i una llavor inicial SEED. Tot i que S té mida 128x128 només les primeres L
+! files i L columnes s'utilitzen.
 subroutine generateSpinMatrix(S, L, SEED)
   use iso_fortran_env
   implicit none
   integer*4, intent(in) :: L
   integer*4, intent(in) :: SEED
-  integer*2, intent(out) :: S(1:64,1:64)
+  integer*2, intent(out) :: S(1:128,1:128)
   integer*4 :: i, j
   integer(int32) :: genrand_int31
 
diff --git a/quad10/fenomens/lab/p4/graphs/dependenciaEnL.gnu b/quad10/fenomens/lab/p4/graphs/dependenciaEnL.gnu
index a1a1ae8..8c04f2c 100755
--- a/quad10/fenomens/lab/p4/graphs/dependenciaEnL.gnu
+++ b/quad10/fenomens/lab/p4/graphs/dependenciaEnL.gnu
@@ -1,104 +1,352 @@
 #!/usr/bin/env -S gnuplot -c
-outputfile = 'data_out/dep_en_L_graphs/' # Nom de la imatge resultant (sense extensió)
-datafilepre = 'data_out/dep_en_L/tmpdata/'
+outputfile = ''
+datafilepre = '../dep_en_L/tmpdata/'
 datafilepost = '.dat'
 LS = ARG1
 
-svgTerminal = "set terminal svg dashed size 800, 800 font \"Computer Modern,Tinos,Helvetica,15\"; pointSize=0.5"
+SPECIAL_L = 32
+ZOOMEDIN = 1 # Whether we want to include zoomed in graphs (1) or not (0)
+USEHACK = 1 # In theory, due to the defaults of run.bash, SUI = 1900*150. But for L=96, we used other values, so setting USEHACK=1 gives a different value for SUMI(96).
+
+# COMPORTAMENT CRÍTIC (extret de l'anàlisi fet a l'informe):
+TC = 2.268
+BETA = 0.13
+GAMMA = 1.77
+
+# All LS except for 96, since it's low quality data.
+PartialLS = system('LS="'. LS .'"; PartialLS=""; for L in $LS; do if [ $L -ne 96 ]; then PartialLS+="$L "; fi; done; echo $PartialLS')
+
+TERMINALS = "svg png tex"
+
+svgTerminal = "set terminal svg size 800, 800 font \"Computer Modern,Tinos,Helvetica,15\"; pointSize=0.5"
 pngTerminal = "set terminal pngcairo size 800, 800 font \"Computer Modern,Tinos,Helvetica,15\"; pointSize=0.7"
+texTerminal = "set terminal cairolatex size 8.5cm, 8.5cm; pointSize=0.7"
 graphStyle = "with linespoints pointsize pointSize"
 
+SUMI(L) = (USEHACK ? (L == 96 ? 1400*16 : 1900*150) : 1900*150)
+suffix(CHOOSEL) = (CHOOSEL == 0 ? '_ALLL' : (CHOOSEL == 1 ? '' : '_L' . SPECIAL_L))
+
 set style line 101 lc rgb '#808080' lt 1 lw 1
 set border 3 front ls 101
 set tics nomirror out scale 0.75
-set style line 102 lc rgb '#d6d7d9' lt 0 lw 1
+set linetype 99963 dashtype (6, 0, 0, 10)
+set style line 102 lc rgb '#d6d7d9' lt 99963 lw 1
 set grid back ls 102
 
-set xlabel "Temperatura (K)"
+set xlabel "T"
 set yrange [*<0:0<*]
 # ------------------------------------
 # Gràfica per la capacitat calorífica:
 # ------------------------------------
-set title "Capacitat calorífica per diferents valors d'L"
+set title ""
 
-@svgTerminal
-set output outputfile.'capacitat_calorifica.svg'
-plot for [L in LS] datafilepre . L . datafilepost using 2:10 @graphStyle title "L = ".L
+# CHOOSEL can be:
+#   - 0: all Ls will be drawn
+#   - 1: all Ls except for 96 will be drawn
+#   - 2: only L=32 will be drawn
+do for [CHOOSEL = 0:2] {
+  do for [TERM in TERMINALS] {
+    if (TERM eq "tex") {
+      @texTerminal
+      set ylabel '$c_V$'
+      title_1 = '$\frac{\text{Var}(e)}{T^2}$'
+      title_1alt = 'L = %d'
+      title_2 = '$\frac{d\langle e \rangle}{dT}$'
+      set key spacing 2
+    } else {
+      set ylabel 'Capacitat calorífica'
+      title_1 = 'Var(e)/T^2'
+      title_1alt = 'L = %d'
+      title_2 = 'd<e>/dT'
+      set key spacing 1.25
+      if (TERM eq "png") {
+        @pngTerminal
+      } else {
+        @svgTerminal
+      }
+    }
 
-@pngTerminal
-set output outputfile.'capacitat_calorifica.png'
-replot
+    set output outputfile.'capacitat_calorifica'.suffix(CHOOSEL).'.'.TERM
+    set yrange [*<0:0<*]
+
+    if (CHOOSEL == 0) {
+      set yrange [0:2.6]
+      plot for [L in LS] datafilepre . L . datafilepost using 2:10 @graphStyle title sprintf(title_1alt, L+0)
+    } else {
+      if (CHOOSEL == 1) {
+        plot for [L in PartialLS] datafilepre . L . datafilepost using 2:10 @graphStyle title sprintf(title_1alt, L+0)
+      } else {
+        plot datafilepre . SPECIAL_L . datafilepost using 2:10 @graphStyle title sprintf(title_1, SPECIAL_L), \
+             datafilepre . SPECIAL_L . '_derivada_E' .datafilepost using 1:($2/(SPECIAL_L**2)) with points pointsize 1.2*pointSize title sprintf(title_2, SPECIAL_L)
+      }
+    }
+  }
+}
+
+set yrange [*<0:0<*]
+set ylabel ""
 
 # ----------------------------------------
 # Gràfica per la susceptibilitat magnètica
 # ----------------------------------------
-set title "Susceptibilitat magnètica per diferents valors d'L"
+set title ""
+
+# CHOOSEL can be:
+#   - 0: all Ls will be drawn
+#   - 1: all Ls except for 96 will be drawn
+#   - 2: only L=32 will be drawn
+do for [CHOOSEL = 0:3] {
+  set ylabel 'Susceptibilitat magnètica'
+
+  @svgTerminal
+  set output outputfile.'susceptibilitat_magnetica'. suffix(CHOOSEL) .'.svg'
+  if (CHOOSEL == 0) {
+    set yrange [0:175]
+    plot for [L in LS] datafilepre . L . datafilepost using 2:11 @graphStyle title "L = ".L
+  } else {
+    if (CHOOSEL == 1) {
+      set yrange [0:75]
+      plot for [L in PartialLS] datafilepre . L . datafilepost using 2:11 @graphStyle title "L = ".L
+    } else {
+      set yrange [0:22]
+      plot datafilepre . SPECIAL_L . datafilepost using 2:11 @graphStyle lc 3 title '$\chi$'
+    }
+  }
+
+  @pngTerminal
+  set output outputfile.'susceptibilitat_magnetica'. suffix(CHOOSEL) .'.png'
+  replot
+
+  @texTerminal
+  set ylabel '$\chi$'
+  set output outputfile.'susceptibilitat_magnetica'. suffix(CHOOSEL) .'.tex'
+  replot
+}
+
+set yrange [*<0:0<*]
+
+# --------------------------------------------------------------
+# Gràfica per la susceptibilitat magnètica (finite size scaling)
+# --------------------------------------------------------------
+set title ""
+
+set xlabel 'L^{1/\nu} t'
+set ylabel '\chi/(L^{\gamma/\nu})'
+
+set key spacing 1.25
 
 @svgTerminal
-set output outputfile.'susceptibilitat_magnetica.svg'
-plot for [L in LS] datafilepre . L . datafilepost using 2:11 @graphStyle title "L = ".L
+set output outputfile.'susceptibilitat_magnetica_fss.svg'
+plot for [L in PartialLS] datafilepre . L . datafilepost using (L*($2 - TC)/TC):($11/(L**GAMMA)) @graphStyle title "L = ".L
 
 @pngTerminal
-set output outputfile.'susceptibilitat_magnetica.png'
+set output outputfile.'susceptibilitat_magnetica_fss.png'
 replot
 
+@texTerminal
+set xlabel '$L^{1/\nu} t$'
+set ylabel '$\chi/L^{\gamma/\nu}$'
+set key spacing 1.66
+set output outputfile.'susceptibilitat_magnetica_fss.tex'
+replot
+
+set key spacing 1.25
+set xlabel 'T'
+
 # ---------------------
 # Gràfica per l'energia
 # ---------------------
-do for [IDX = 0:1] {
-  if (IDX == 0) {
-    @svgTerminal
-    set output outputfile.'energia.svg'
-  } else {
-    @pngTerminal
-    set output outputfile.'energia.png'
+
+# CHOOSEL can be:
+#   - 0: all Ls will be drawn
+#   - 1: all Ls except for 96 will be drawn
+#   - 2: only L=32 will be drawn
+do for [CHOOSEL = 0:2] {
+  do for [TERM in TERMINALS] {
+    if (TERM eq "tex") {
+      @texTerminal
+      title_l = 'L = %d'
+      title_e = '$\langle e \rangle$'
+      title_sqrtee = '$-\sqrt{\langle e^2 \rangle}$'
+      set key spacing 1.66
+    } else {
+      title_l = 'L = %d'
+      title_e = '<e>'
+      title_sqrtee = '-sqrt(<e^2>)'
+      set key spacing 1.25
+      if (TERM eq "png") {
+        @pngTerminal
+      } else {
+        @svgTerminal
+      }
+    }
+    set output outputfile.'energia'.suffix(CHOOSEL).'.'.TERM
+
+    set origin 0, 0
+    set size 1, 1
+    set xtics auto
+    set ytics auto
+
+    if (ZOOMEDIN = 1) {
+      set multiplot
+    }
+
+    unset object 1
+
+    set title ""
+    set xlabel "T"
+    if (TERM eq "tex" && CHOOSEL != 2) {
+      set ylabel '$\langle e \rangle$'
+    } else {
+      set ylabel "Energia per partícula"
+    }
+    set key bottom right
+    set key noopaque
+    set autoscale x
+    set yrange [*<0:0<*]
+
+    if (TERM eq "tex") {
+      set ytics 0.5
+    }
+
+    if (CHOOSEL == 0) {
+      plot for [L in LS] datafilepre . L . datafilepost using 2:($3/(L**2)) @graphStyle title sprintf(title_l, L+0)
+    } else {
+      if (CHOOSEL == 1) {
+        plot for [L in PartialLS] datafilepre . L . datafilepost using 2:($3/(L**2)) @graphStyle title sprintf(title_l, L+0)
+      } else {
+        plot datafilepre . SPECIAL_L . datafilepost using 2:($3/(SPECIAL_L**2)) @graphStyle title sprintf(title_e, SPECIAL_L), \
+             datafilepre . SPECIAL_L . datafilepost using 2:(-sqrt($4)/(SPECIAL_L**2)) @graphStyle title sprintf(title_sqrtee, SPECIAL_L)
+      }
+    }
+
+    if (ZOOMEDIN = 1) {
+      if (TERM eq "tex") {
+        set origin 0.175, 0.575
+      } else {
+        set origin 0.1, 0.525
+      }
+      set size 0.6, 0.4
+      set title ""
+      set xlabel ""
+      set ylabel ""
+      set key off
+      if (CHOOSEL != 2) {
+        set xrange [2.3:2.36]
+        set yrange [-1.5:-1.2]
+        if (TERM eq "tex") {
+          set xtics 0.03
+          set ytics 0.15
+        }
+      } else {
+        set xrange [2.35:2.4]
+        set yrange [-1.3:-1.2]
+        set xtics 0.025
+        if (TERM eq "tex") {
+          set ytics 0.05
+        }
+      }
+      set object 1 rectangle from graph 0,0 to graph 1,1 behind fillcolor rgb 'white' fillstyle solid noborder
+      replot
+
+      unset multiplot
+    }
   }
-
-  set multiplot
-
-  set origin 0, 0
-  set size 1, 1
-  set title "Energia per diferents valors de L"
-  set xlabel "Temperatura (K)"
-  set key bottom right
-  set key noopaque
-  set autoscale x
-  set yrange [*<0:0<*]
-  unset object 1
-  plot for [L in LS] datafilepre . L . datafilepost using 2:($3/(L**2)) @graphStyle title "<E>/N, L = ".L, \
-       for [L in LS] datafilepre . L . datafilepost using 2:(-sqrt($4)/(L**2)) @graphStyle title "-sqrt(<E^2>)/N, L = ".L
-
-  set origin 0.1, 0.525
-  set size 0.6, 0.4
-  set title ""
-  set xlabel ""
-  set key off
-  set xrange [2.25:2.5]
-  set yrange [-1.6:-1.1]
-  set object 1 rectangle from graph 0,0 to graph 1,1 behind fillcolor rgb 'white' fillstyle solid noborder
-  replot
-
-  unset multiplot
 }
 
+set key spacing 1.25
 set origin 0, 0
 set size 1, 1
 set autoscale x
 set yrange [*<0:0<*]
-set xlabel "Temperatura (K)"
+set xlabel "T"
+set ylabel ""
+set xtics auto
+set ytics auto
 unset object 1
 
 # ----------------------------
 # Gràfica per la magnetització
 # ----------------------------
-set title "Magnetització per diferents valors d'L"
+set title ""
 
-@svgTerminal
-set output outputfile.'magnetitzacio.svg'
-set key top right
-plot for [L in LS] datafilepre . L . datafilepost using 2:($7/(L**2)) @graphStyle title "<|M|>/N, L = ".L, \
-     for [L in LS] datafilepre . L . datafilepost using 2:(sqrt($8)/(L**2)) @graphStyle title "sqrt(<M^2>)/N, L = ".L
+# CHOOSEL can be:
+#   - 0: all Ls will be drawn
+#   - 1: all Ls except for 96 will be drawn
+#   - 2: only L=32 will be drawn
+do for [CHOOSEL = 0:2] {
+  do for [TERM in TERMINALS] {
+    if (TERM eq "tex") {
+      @texTerminal
+      title_m = '$\langle |m| \rangle$'
+      title_malt = '$L = %d$'
+      title_sqrtmm = '$\sqrt{\langle m^2 \rangle}$'
+      if (CHOOSEL == 2) {
+        set key spacing 2
+      } else {
+        set key spacing 1.66
+      }
+    } else {
+      title_m = '<|m|>'
+      title_malt = 'L = %d'
+      title_sqrtmm = 'sqrt(<m^2>)'
+      set key spacing 1.25
+      if (TERM eq "png") {
+        @pngTerminal
+      } else {
+        @svgTerminal
+      }
+    }
 
-@pngTerminal
-set output outputfile.'magnetitzacio.png'
-replot
+    if (TERM eq "tex" && CHOOSEL != 2) {
+      set ylabel '$\langle |m| \rangle$'
+    } else {
+      set ylabel "Magnetització per partícula"
+    }
+
+    set output outputfile.'magnetitzacio'.suffix(CHOOSEL).'.'.TERM
+    set key top right
+    if (CHOOSEL == 0) {
+      plot for [L in LS] datafilepre . L . datafilepost using 2:($7/(L**2)) @graphStyle title sprintf(title_malt, L+0)
+    } else {
+      if (CHOOSEL == 1) {
+        plot for [L in PartialLS] datafilepre . L . datafilepost using 2:($7/(L**2)) @graphStyle title sprintf(title_malt, L+0)
+      } else {
+        plot datafilepre . SPECIAL_L . datafilepost using 2:($7/(SPECIAL_L**2)) @graphStyle title sprintf(title_m, SPECIAL_L), \
+             datafilepre . SPECIAL_L . datafilepost using 2:(sqrt($8)/(SPECIAL_L**2)) @graphStyle title sprintf(title_sqrtmm, SPECIAL_L)
+      }
+    }
+  }
+}
+
+# --------------------------------------------------
+# Gràfica per la magnetització (finite size scaling)
+# --------------------------------------------------
+set title ""
+set ylabel "Magnetització per partícula"
+
+# CHOOSEL can be:
+#   - 0: all Ls will be drawn
+#   - 1: all Ls except for 96 will be drawn
+#   - 2: only L=32 will be drawn
+do for [TERM in TERMINALS] {
+  if (TERM eq "tex") {
+    @texTerminal
+    set xlabel '$L^{1/\nu} t$'
+    set ylabel '$\langle |m| \rangle L^{\beta/\nu}$'
+    set key spacing 1.66
+  } else {
+    set xlabel 'L^{1/\nu} t'
+    set ylabel '<|m|> L^{\beta/\nu}'
+    set key spacing 1.25
+    if (TERM eq "png") {
+      @pngTerminal
+    } else {
+      @svgTerminal
+    }
+  }
+
+  set output outputfile.'magnetitzacio_fss.'.TERM
+  set key top right
+  plot for [L in PartialLS] datafilepre . L . datafilepost using (L*($2 - TC)/TC):($7/(L**2)*(L**(BETA))) @graphStyle title sprintf('L = %s', L)
+}
diff --git a/quad10/fenomens/lab/p4/graphs/dependenciaEnSeed.gnu b/quad10/fenomens/lab/p4/graphs/dependenciaEnSeed.gnu
new file mode 100755
index 0000000..be9b69c
--- /dev/null
+++ b/quad10/fenomens/lab/p4/graphs/dependenciaEnSeed.gnu
@@ -0,0 +1,31 @@
+#!/usr/bin/env -S gnuplot -c
+outputfile = 'depEnSeed' # Nom de la imatge resultant (sense extensió)
+datafilepost = '.ev'
+L=32
+SEEDS = ARG1
+L = ARG2
+
+do for [IDX = 0:1] {
+  if (IDX == 0) {
+    set terminal cairolatex size 8.5cm, 11cm
+    set output outputfile.'.tex'
+  } else {
+    set terminal svg dashed size 600, 1200 font "Computer Modern,Tinos,Helvetica,15"
+    set output outputfile.'.svg'
+  }
+
+  set multiplot layout 2,1
+  set key outside top horizontal
+
+  set xlabel "Iteracions MC"
+  set ylabel '$\langle e \rangle$'
+  set yrange [-2:0]
+  plot for [S in SEEDS] S . datafilepost using 1:($2/(L*L)) with lines title sprintf('$S = %s$', S)
+
+  set xlabel "Iteracions MC"
+  set ylabel '$\langle m \rangle$'
+  set yrange [-1:1]
+  plot for [S in SEEDS] S . datafilepost using 1:($3/(L*L)) with lines title sprintf('$S = %s$', S)
+
+  unset multiplot
+}
diff --git a/quad10/fenomens/lab/p4/graphs/finalConfs.gnu b/quad10/fenomens/lab/p4/graphs/finalConfs.gnu
new file mode 100755
index 0000000..c83fa5f
--- /dev/null
+++ b/quad10/fenomens/lab/p4/graphs/finalConfs.gnu
@@ -0,0 +1,31 @@
+#!/usr/bin/env -S gnuplot -c
+outputfile = 'confs' # Nom de la imatge resultant (sense extensió)
+datafilepost = '.conf'
+TS = ARG1
+L = ARG2
+
+do for [IDX = 0:1] {
+  if (IDX == 0) {
+    set terminal cairolatex size 8.5cm, 13cm
+    set output outputfile.'.tex'
+    symbsize=0.18
+  } else {
+    set terminal svg dashed size 600, 1000 font "Computer Modern,Tinos,Helvetica,15"
+    set output outputfile.'.svg'
+    symbsize=0.3
+  }
+
+  set multiplot layout 3,2
+
+  do for [T in TS] {
+    set title "T = " . T
+    set size square
+    set xrange [0:L+1]
+    set yrange [0:L+1]
+    unset xtics
+    unset ytics
+    plot T . datafilepost using 1:2 with points pt 5 ps symbsize notitle
+  }
+
+  unset multiplot
+}
diff --git a/quad10/fenomens/lab/p4/magne.f90 b/quad10/fenomens/lab/p4/magne.f90
index 118a965..4e0e89e 100644
--- a/quad10/fenomens/lab/p4/magne.f90
+++ b/quad10/fenomens/lab/p4/magne.f90
@@ -1,7 +1,9 @@
+! Retorna, donada una matriu d'spins S d'un sistema de mida LxL,
+! la magnetització del sistema.
 real*8 function magne(S, L)
   implicit none
   integer*4, intent(in) :: L
-  integer*2, intent(out) :: S(1:64,1:64)
+  integer*2, intent(out) :: S(1:128,1:128)
   integer*4 :: i, j
   real*8 :: mag
   mag = 0d0
diff --git a/quad10/fenomens/lab/p4/pbc.f90 b/quad10/fenomens/lab/p4/pbc.f90
index 09cfff3..3f803a4 100644
--- a/quad10/fenomens/lab/p4/pbc.f90
+++ b/quad10/fenomens/lab/p4/pbc.f90
@@ -1,3 +1,7 @@
+! Funció que imposa les condicions de contorn periòdiques, enviant l'índex N al
+! valor corresponent en l'interval [1, L].
+!
+! Prerequisit: N està a l'interval [0, L + 1].
 integer*4 function PBC(N, L)
   integer*4, intent(in) :: N, L
   if (N == 0) then
diff --git a/quad10/fenomens/lab/p4/run.bash b/quad10/fenomens/lab/p4/run.bash
index 2c95bdd..1be282a 100755
--- a/quad10/fenomens/lab/p4/run.bash
+++ b/quad10/fenomens/lab/p4/run.bash
@@ -5,9 +5,10 @@
 
   Usage: $progname T L [outFilePrefix [--help --skipIfComputed
              --nSeeds NSEEDS --mcTot MCTOT --mcIni MCINI --mcD MCD
-             --initialSeed INITIALSEED --saveFinalConf]]
+             --initialSeed INITIALSEED --saveFinalConf
+             --saveEvolution]]
 
-  the output file will be saved at "data_out/{{outFilePrefix}}"
+  the results file will be saved at "data_out/{{outFilePrefix}}.res"
 
   optional arguments:
     -h, --help            show this help message and exit.
@@ -24,6 +25,11 @@
     --initialSeed         initial seed for the Monte Carlo algorithm.
     --saveFinalConf       save a file with the configuration which
                           belongs to the last Montecarlo iteration.
+                          It will be saved at
+                          "data_out/{{outFilePrefix}}.conf".
+    --saveEvolution       save a file with the temporary evolution
+                          of the simulation. It will be saved at
+                          "data_out/{{outFilePrefix}}.ev"
 END
 }
 
@@ -37,7 +43,7 @@
 defaultOutFilePrefix="L$L-T$T"
 outFilePrefix="${3:-$defaultOutFilePrefix}"
 
-opts=$(getopt -l "help,skipIfComputed,nSeeds:,mcTot:,mcIni:,mcD:,initialSeed:,saveFinalConf" -o "hsn:m:i:d:" -n "$progname" -- "${@:4}")
+opts=$(getopt -l "help,skipIfComputed,nSeeds:,mcTot:,mcIni:,mcD:,initialSeed:,saveFinalConf,saveEvolution" -o "hsn:m:i:d:" -n "$progname" -- "${@:4}")
 eval set -- "$opts"
 
 skipIfComputed=false
@@ -47,6 +53,7 @@
 mcD=20
 initialSeed=117654
 saveFinalConf=F
+saveEvolution=F
 
 while true; do
   case "$1" in
@@ -82,6 +89,10 @@
       saveFinalConf=T
       shift
       ;;
+    --saveEvolution)
+      saveEvolution=T
+      shift
+      ;;
     *) break ;;
   esac
 done
@@ -104,5 +115,6 @@
 MCINI=$mcIni,
 MCD=$mcD
 SAVEFINALCONF=$saveFinalConf
+SAVEEVOLUTION=$saveEvolution
 &END
 EOF
diff --git a/quad10/fenomens/lab/p4/runAll.bash b/quad10/fenomens/lab/p4/runAll.bash
index c4f0783..2a842a0 100755
--- a/quad10/fenomens/lab/p4/runAll.bash
+++ b/quad10/fenomens/lab/p4/runAll.bash
@@ -1,4 +1,4 @@
-TList=$((seq 1.6 0.05 3.2; seq 2.2 0.01 2.5) | sort -n | uniq)
+TList=$((seq 1.6 0.05 3.2; seq 2.2 0.01 2.57; seq 2.25 0.005 2.4; seq 2.39 0.005 2.41; seq 2.53 0.005 2.55) | sort -n | uniq)
 LList="8 16 32 64"
 for L in $LList; do
   mkdir -p data_out/dep_en_L/$L
diff --git a/quad10/fenomens/lab/p4/writeconfig.f90 b/quad10/fenomens/lab/p4/writeconfig.f90
index 30de3c9..477129f 100644
--- a/quad10/fenomens/lab/p4/writeconfig.f90
+++ b/quad10/fenomens/lab/p4/writeconfig.f90
@@ -1,3 +1,5 @@
+! Funció que guarda la configuració S d'un sistema de mida LxL al fitxer amb
+! nom NOM.
 subroutine writeConfig(S, L, NOM)
   implicit none
   integer*4, intent(in) :: L