diff --git a/functions.py b/functions.py
index 28265d5..0426034 100644
--- a/functions.py
+++ b/functions.py
@@ -12,7 +12,7 @@ from scipy.signal import find_peaks
 import matplotlib.transforms as mtransforms
 
 
-def error_plot(folder,t_step,r_criteria,save):
+def error_plot(folder,t_step,r_criteria,save_path):
     # Load iterations and residual error
     histor = folder + '/histor.dat'
     input = open(histor, 'r')
@@ -29,30 +29,31 @@ def error_plot(folder,t_step,r_criteria,save):
             error.append(error_info[2][n-1])
 
     time = np.linspace(start = t_step, stop = len(error)*t_step, num = len(error))
-    
-    # Plots of interest
+
     # Liniar Scale
-    plt.figure()  
-    plt.plot(time,error)
-    plt.plot(time,r_criteria*np.ones(len(error)),'r')
-    plt.ylabel('Residual error')
-    plt.xlabel('Time steps')
-    plt.title('Last nonlinear residual error for each time step')
-    plt.grid(True)
-    if save: plt.savefig(plt_folder + case + '_Last_nonlin_res_error.pdf')
-
+    fig, ax = plt.subplots()
+    ax.plot(time,error)
+    ax.plot(time,r_criteria*np.ones(len(error)),'r')
+    ax.set(xlabel='Time steps', ylabel='Residual error',
+        title='Last nonlinear residual error for each time step')
+    ax.spines['right'].set_visible(False)
+    ax.spines['top'].set_visible(False)
+    plt.show()
+    plt.savefig(save_path + '/Last_nonlin_res_error.pdf')
     # Semilog scale
-    plt.figure()  
-    plt.semilogy(time,error)
-    plt.semilogy(time,r_criteria*np.ones(len(error)),'r')
-    plt.ylabel('Residual error')
-    plt.xlabel('Time steps')
-    plt.title('Log - Last nonlinear residual error for each time step')
-    plt.grid(True)
-    if save: plt.savefig(plt_folder + case + '_Log_Last_nonlin_res_error.pdf')
-    
+    fig, ax = plt.subplots()
+    ax.semilogy(time,error)
+    ax.semilogy(time,r_criteria*np.ones(len(error)),'r')
+    ax.set(xlabel='Time steps', ylabel='Residual error',
+        title='Log - Last nonlinear residual error for each time step')
+    ax.spines['right'].set_visible(False)
+    ax.spines['top'].set_visible(False)
+    plt.show()
+    plt.savefig(save_path + '/Log_Last_nonlin_res_error.pdf')
+    fig.savefig(save_path + '/Log_Last_nonlin_res_error.jpg',dpi=150)
 
-def periodicity(project,folder,dt,T_cyc,n_cyc):
+
+def periodicity(project,folder,dt,T_cyc,n_cyc,save_path):
     pressure = np.loadtxt(folder+'/PHistRCR.dat',skiprows=2,)
     time = np.linspace(0,T_cyc*n_cyc,round(T_cyc/dt*n_cyc))
     peak_P = []
@@ -66,19 +67,23 @@ def periodicity(project,folder,dt,T_cyc,n_cyc):
     peak_Pdiff = list(map(abs, peak_Pdiff)) 
     
     fig, ax = plt.subplots()
-    ax.plot(time,pressure[:,-1]/1333.22)
+    ax.plot(time,pressure[:,-1]/1333.22,'b')
     ax.plot(time[peak_P_pos], peak_P,'ro',label='Cylce pike')
     ax.set(xlabel='time [s]', ylabel='Pressure [mmHg]',
-        title='Pressure @ last outlet')
+        title='Pressure at last outlet')
     ax.spines['right'].set_visible(False)
     ax.spines['top'].set_visible(False)
     ax.legend(loc=0)
     plt.show()
+    plt.savefig(save_path + '/periodicity.pdf')
+    fig.savefig(save_path + '/periodicity.jpg',dpi=150)
     
-    if (peak_Pdiff[-1]<=1): print('The numerical simulation \'{0}\' has achieve periodicity!\nSystolic Blood Pressure (SBP):\nsecond-last cycle = {1:.2f} mmHg,\nlast cycle = {2:.2f} mmHg,\n\u0394mmHg = {3:.2f} mmHg'.format(project,peak_P[-2],peak_P[-1],peak_Pdiff[-1]))  
-
-
-def pressure(folder,N_ts,T_cyc,dt,n_cyc):
+    if (peak_Pdiff[-1]<=1): 
+        print('The numerical simulation \'{0}\' has achieve periodicity!\nSystolic Blood Pressure (SBP):\nsecond-last cycle = {1:.2f} mmHg,\nlast cycle = {2:.2f} mmHg,\n\u0394mmHg = {3:.2f} mmHg'.format(project,peak_P[-2],peak_P[-1],peak_Pdiff[-1]))  
+        txt = 'The numerical simulation \'{0}\' has achieve periodicity!\nSystolic Blood Pressure (SBP):\nsecond-last cycle = {1:.2f} mmHg,\nlast cycle = {2:.2f} mmHg,\nDelta_mmHg = {3:.2f} mmHg'.format(project,peak_P[-2],peak_P[-1],peak_Pdiff[-1])
+    return txt
+    
+def pressure(folder,N_ts,T_cyc,dt,n_cyc,save_path):
     pressure = np.loadtxt(folder+'/PHistRCR.dat',skiprows=2,)
     Nc = round(T_cyc/dt)
     time = np.linspace(0,T_cyc,Nc)
@@ -93,15 +98,17 @@ def pressure(folder,N_ts,T_cyc,dt,n_cyc):
         MBP[i] = (mean(pressure[N_ts-Nc:N_ts,i]/1333.22)) 
     PP = SBP-DBP    
     ax.set(xlabel='time [s]', ylabel='Pressure [mmHg]',
-        title='Pressure @ each outlet')
+        title='Pressure at each outlet')
     ax.spines['right'].set_visible(False)
     ax.spines['top'].set_visible(False)
     ax.legend(loc=0)
     plt.show()
+    plt.savefig(save_path + '/pressure.pdf')
+    fig.savefig(save_path + '/pressure.jpg',dpi=150)
     return (DBP,MBP,SBP,PP)
 
 
-def flow(folder,N_ts,T_cyc,dt,n_cyc):
+def flow(folder,N_ts,T_cyc,dt,n_cyc,save_path):
     flow = np.loadtxt(folder+'/QHistRCR.dat',skiprows=2,)
     Nc = round(T_cyc/dt)
     time = np.linspace(0,T_cyc,Nc)
@@ -112,14 +119,16 @@ def flow(folder,N_ts,T_cyc,dt,n_cyc):
         Q[i] = (mean(flow[N_ts-Nc:N_ts,i]))
    
     ax.set(xlabel='time [s]', ylabel='Flow [mL/s]',
-        title='Flow @ each outlet')
+        title='Flow at each outlet')
     ax.spines['right'].set_visible(False)
     ax.spines['top'].set_visible(False)
     ax.legend(loc=0)
     plt.show()
+    plt.savefig(save_path + '/flow.pdf')
+    fig.savefig(save_path + '/flow.jpg',dpi=150)
     return Q
     
-def inlet_flow_waveform(project_folder,t_btw_rst,N_ts,dt,T_cyc,n_cyc):
+def inlet_flow_waveform(project_folder,t_btw_rst,N_ts,dt,T_cyc,n_cyc,save_path):
     x = np.loadtxt(project_folder+'/ROI-1.flow')
     t = x[:,0]
     Q = -x[:,1]
@@ -136,12 +145,11 @@ def inlet_flow_waveform(project_folder,t_btw_rst,N_ts,dt,T_cyc,n_cyc):
 
     fig, ax = plt.subplots()
     ax.plot(t, Q, 'r')
-    ax.plot(t_pts, Q_pts, 'ob')
+    ax.plot(t_pts, Q_pts, 'ob',label='Time steps saved')
     trans_offset = mtransforms.offset_copy(ax.transData, fig=fig,
                                        x=-0, y=0.15, units='inches')
-
     ax.set(xlabel='Time [s]', ylabel='Flow Rate - Q [mL/s]',
-       title='Inlet Flow rate Waveform - 1 cycle')
+       title='Inlet Flow Rate Waveform')
     ax.set_ylim([-10, 90])
     ax.spines['right'].set_visible(False)
     ax.spines['top'].set_visible(False)
@@ -151,8 +159,14 @@ def inlet_flow_waveform(project_folder,t_btw_rst,N_ts,dt,T_cyc,n_cyc):
     for i in range(0,np.unique(np.round(t_pts,3)).shape[0]):
         time.append('$t_'+str(i+1)+'$')
     for x, y, t in zip(t_pts[(-n_cyc-1):], Q_pts[(-n_cyc-1):], time):
-        plt.text(x, y, t, transform=trans_offset, fontsize=12)        
+        plt.text(x, y, t, transform=trans_offset, fontsize=12)    
+    ax.legend(loc=0)
     plt.show()
-            
+    plt.savefig(save_path + '/inlet_waveform.pdf')
+    fig.savefig(save_path + '/inlet_waveform.jpg',dpi=150)
+    print('{0} time steps saved, available to visualize in ParaView.'.format(np.unique(np.round(t_pts,3)).shape[0]))
+    txt = '{0} time steps saved, available to visualize in ParaView.'.format(np.unique(np.round(t_pts,3)).shape[0])
+    return txt
+      
         
         
\ No newline at end of file
diff --git a/main.py b/main.py
index 0de3baa..e818514 100644
--- a/main.py
+++ b/main.py
@@ -11,6 +11,7 @@ import os.path
 from scipy.signal import find_peaks_cwt
 from scipy import signal
 import statistics
+from fpdf import FPDF
 
 
 from functions import error_plot, periodicity, pressure, flow, inlet_flow_waveform
@@ -20,6 +21,9 @@ Tk().withdraw()
 folder = askdirectory()
 project_folder = os.path.dirname(folder)
 project = os.path.basename(project_folder)
+save_path = folder+'/'+project+'-report'
+# os.mkdir(save_path)
+save_pdf = save_path + '/' + project + '-report.pdf'
 
 # Input parameters
 T_cyc = 0.477
@@ -42,14 +46,36 @@ rc = float(mylines[26][18:-1])
 t_btw_rst = int(mylines[6][37:-1])
 
 # Cehcking convergency and periodicity
-error_plot(folder,dt,rc,False)
-periodicity(project,folder,dt,T_cyc,n_cyc)
+error_plot(folder,dt,rc,save_path)
+txt1 = periodicity(project,folder,dt,T_cyc,n_cyc,save_path)
 
 # Pressure - Outlets
-(DBP,MBP,SBP,PP) = pressure(folder,N_ts,T_cyc,dt,n_cyc)
+(DBP,MBP,SBP,PP) = pressure(folder,N_ts,T_cyc,dt,n_cyc,save_path)
  
 # Flow Rate - Outlets
-(Q_avg) = flow(folder,N_ts,T_cyc,dt,n_cyc)
+(Q_avg) = flow(folder,N_ts,T_cyc,dt,n_cyc,save_path)
 
 # Inlet Flow Rate + and t saved
-inlet_flow_waveform(project_folder,t_btw_rst,N_ts,dt,T_cyc,n_cyc)
\ No newline at end of file
+txt2 = inlet_flow_waveform(project_folder,t_btw_rst,N_ts,dt,T_cyc,n_cyc,save_path)
+
+
+# Create PDF report
+pdf = FPDF('P','mm','Letter')
+pdf.add_page()
+pdf.set_font('Times', 'B', 18)
+pdf.cell(200, 20, 'Project name: '+ project, 0, 1, 'C')
+
+# Keep this ratio - 5.42/3.86
+pdf.set_font('Times', '', 14)
+pdf.cell(200, 30, 'Cehcking convergency and periodicity:', 0, 1, 'L')
+pdf.image(save_path +'/Log_Last_nonlin_res_error.jpg', x = None, y = None, w = 140, h = 100, type = '', link = '')
+pdf.image(save_path +'/periodicity.jpg', x = None, y = None, w = 140, h = 100, type = '', link = '')
+pdf.cell(200, 20, 'Cehcking Pressures at each outlet:', 0, 1, 'L')
+pdf.image(save_path +'/pressure.jpg', x = None, y = None, w = 140, h = 100, type = '', link = '')
+pdf.cell(200, 20, 'Cehcking Flow Rate at each outlet:', 0, 1, 'L')
+pdf.image(save_path +'/flow.jpg', x = None, y = None, w = 140, h = 100, type = '', link = '')
+
+pdf.image(save_path +'/inlet_waveform.jpg', x = None, y = None, w = 140, h = 100, type = '', link = '')
+
+
+pdf.output(save_pdf, 'F')