diff --git a/wide vs narrow/PT curve/draw_pt_curve.m b/wide vs narrow/PT curve/draw_pt_curve.m
new file mode 100755
index 0000000..da896b9
--- /dev/null
+++ b/wide vs narrow/PT curve/draw_pt_curve.m	
@@ -0,0 +1,11 @@
+clc; clear; close all;
+
+N = 256;
+M = 4;
+
+[K_range, N_b, N_s] = get_phase_transition_curve(N, M);
+
+figure;
+plot(K_range, N_b);
+hold on;
+plot(K_range, N_b);
diff --git a/wide vs narrow/PT curve/get_integral.m b/wide vs narrow/PT curve/get_integral.m
new file mode 100755
index 0000000..dd2a3ae
--- /dev/null
+++ b/wide vs narrow/PT curve/get_integral.m	
@@ -0,0 +1,20 @@
+function I = get_integral(cache_filename, var_name, tau_range, M)
+    syms x f;
+
+    if exist(cache_filename, "file")
+        data = load(cache_filename, var_name);
+        I = data.(var_name);
+    else
+        I = zeros(length(tau_range), 0);
+
+        for tau_idx = 1:length(tau_range)
+            tau = tau_range(tau_idx);
+            f = (x - tau) ^ 2 * exp(-x ^ 2/2) * x ^ (M - 1) / (2 ^ (M / 2 - 1) * gamma(M / 2));
+            I(tau_idx) = double(int(f, [tau, +inf]));
+        end
+
+        S.(var_name) = I;
+        save(cache_filename, '-struct', 'S', var_name);
+    end
+
+end
diff --git a/wide vs narrow/PT curve/get_phase_transition_curve.m b/wide vs narrow/PT curve/get_phase_transition_curve.m
new file mode 100755
index 0000000..50f39b3
--- /dev/null
+++ b/wide vs narrow/PT curve/get_phase_transition_curve.m	
@@ -0,0 +1,47 @@
+function [K_range, N_b, N_s] = get_phase_transition_curve(N, M, K_max)
+    % This function give the phase transition curve of FAR
+    % input: N, M: FAR parameters
+    % output: N_b, N_s: the phase transition curve of FAR in block case and scalar case
+
+    if nargin < 3
+        K_max = 25;
+    end
+
+    tau_min = 0;
+    tau_max = 15;
+    tau_interval = 0.05;
+    tau_range = tau_min:tau_interval:tau_max;
+
+    K_range = 1:K_max;
+
+    N_b = zeros(length(K_range), 1);
+    N_s = zeros(length(K_range), 1);
+
+    cache_filename_1 = ...
+        "../data/PT_curve_data/" + ...
+        "I_2.mat";
+
+    cache_filename_2 = ...
+        "../data/PT_curve_data/" + ...
+        "I_" + "" + string(2 * M) + ...
+        ".mat";
+
+    I_2 = get_integral(cache_filename_1, "I_2", tau_range, 2);
+    I_2M = get_integral(cache_filename_2, "I_2M", tau_range, 2 * M);
+
+    for K_idx = 1:length(K_range)
+        K = K_range(K_idx);
+        f_set_1 = zeros(length(tau_range), 1);
+        f_set_2 = zeros(length(tau_range), 1);
+
+        for tau_idx = 1:length(tau_range)
+            tau = tau_range(tau_idx);
+            f_set_1(tau_idx) = 1/2 * (K * (2 * M + tau ^ 2) + (N - K) * I_2M(tau_idx));
+            f_set_2(tau_idx) = M / 2 * (K * (2 + tau ^ 2) + (N - K) * I_2(tau_idx));
+        end
+
+        N_b(K_idx) = min(f_set_1);
+        N_s(K_idx) = min(f_set_2);
+    end
+
+end