diff --git a/docs/source/algorithms/CropWorkspace-v1.rst b/docs/source/algorithms/CropWorkspace-v1.rst
index 0257d69e956dd04e0dd827eadb24ff957893a34c..7723bed4a499edeb4fcb81e1da4f29f78b73ec75 100644
--- a/docs/source/algorithms/CropWorkspace-v1.rst
+++ b/docs/source/algorithms/CropWorkspace-v1.rst
@@ -39,8 +39,8 @@ Usage
    OutputWorkspace = CropWorkspace(InputWorkspace=ws,XMin=10.0,XMax=40.0)
 
    # Show workspaces
-   print "TOF Before CropWorkspace",ws.readX(0)
-   print "TOF After CropWorkspace",OutputWorkspace.readX(0)
+   print("TOF Before CropWorkspace {}".format(ws.readX(0)))
+   print("TOF After CropWorkspace {}".format(OutputWorkspace.readX(0)))
    
 Output:
 
diff --git a/docs/source/algorithms/CrossCorrelate-v1.rst b/docs/source/algorithms/CrossCorrelate-v1.rst
index c906e32f125cf955f92c6a9f7279a5701a0fa3f8..cfd7452b8dc7aa1830ab09ee8817daa96d7b23eb 100644
--- a/docs/source/algorithms/CrossCorrelate-v1.rst
+++ b/docs/source/algorithms/CrossCorrelate-v1.rst
@@ -35,8 +35,8 @@ Usage
    OutputWorkspace = CrossCorrelate(InputWorkspace='ws', WorkspaceIndexMax=1, XMin=2, XMax=4)
 
    # Show workspaces
-   print "AutoCorrelation",OutputWorkspace.readY(0)
-   print "CrossCorrelation",OutputWorkspace.readY(1)
+   print("AutoCorrelation {}".format(OutputWorkspace.readY(0)))
+   print("CrossCorrelation {}".format(OutputWorkspace.readY(1)))
 
 .. testoutput:: ExCrossCorrelate
 
diff --git a/docs/source/algorithms/CrystalFieldEnergies-v1.rst b/docs/source/algorithms/CrystalFieldEnergies-v1.rst
index d72c6ca1fc3baf713bdc040a6348f6b851c5d442..cf102a0f030c99b473faf6398231c0a3658f7337 100644
--- a/docs/source/algorithms/CrystalFieldEnergies-v1.rst
+++ b/docs/source/algorithms/CrystalFieldEnergies-v1.rst
@@ -21,11 +21,11 @@ The algorithm calculates the crystal field energies and wave functions. The exam
       en, wf, ham = energies(1,  B20=0.37737, B22=3.9770, B40=-0.031787, B42=-0.11611, B44=-0.12544)
       
       # a list of crystal field energies
-      print 'energies:\n', en
+      print('energies:\n{}'.format(en))
       # a complex-valued matrix with wave functions
-      print 'wave functions:\n', wf
+      print('wave functions:\n{}'.format(wf))
       # a complex-valued matrix with the Hamiltonian
-      print 'Hamiltonian:\n', ham
+      print('Hamiltonian:\n{}'.format(ham))
 
 .. testoutput::
 
diff --git a/docs/source/algorithms/CuboidGaugeVolumeAbsorption-v1.rst b/docs/source/algorithms/CuboidGaugeVolumeAbsorption-v1.rst
index a5d91184e46a4493f8d207e3dc622d8dba08966d..8b5b26683731a737ad4f26c1d4a02ada33cd63f2 100644
--- a/docs/source/algorithms/CuboidGaugeVolumeAbsorption-v1.rst
+++ b/docs/source/algorithms/CuboidGaugeVolumeAbsorption-v1.rst
@@ -65,7 +65,7 @@ Usage
     wsOut = CuboidGaugeVolumeAbsorption(ws, NumberOfWavelengthPoints=5, ElementSize=3,
         SampleHeight=1,SampleWidth=2,SampleThickness=3)
 
-    print "The created workspace has one entry for each spectra: %i" % wsOut.getNumberHistograms()
+    print("The created workspace has one entry for each spectra: {}".format(wsOut.getNumberHistograms()))
 
 Output:
 
diff --git a/docs/source/algorithms/CutMD-v1.rst b/docs/source/algorithms/CutMD-v1.rst
index 6da38bdbf8aedb594ddc2a8c487bd41550b7eefa..7c77511cbc1199e97c023231229cbcc71867c1f1 100644
--- a/docs/source/algorithms/CutMD-v1.rst
+++ b/docs/source/algorithms/CutMD-v1.rst
@@ -157,11 +157,11 @@ _`Usage`
    #Another way we can call CutMD:
    #[out1, out2, out3] = CutMD([to_cut, "some_other_file.nxs", "some_workspace_name"], ...)
 
-   print 'number of dimensions', out_md.getNumDims()
-   print 'number of dimensions not integrated', len(out_md.getNonIntegratedDimensions())
+   print('number of dimensions {}'.format(out_md.getNumDims()))
+   print('number of dimensions not integrated {}'.format(len(out_md.getNonIntegratedDimensions())))
    dim_dE = out_md.getDimension(3)
-   print 'min dE', dim_dE.getMaximum()
-   print 'max dE', dim_dE.getMinimum()
+   print('min dE {}'.format(dim_dE.getMaximum()))
+   print('max dE {}'.format(dim_dE.getMinimum()))
 
 Output:
 
@@ -184,15 +184,15 @@ Output:
    # Cut the MDHistoWorkspace to give a single bin containing half the data              
    cut= CutMD(InputWorkspace=histo_ws, PBins=[[-10, 10], [-5, 5]]) 
 
-   print 'Total signal in input = %0.2f' %  sum(signal)
-   print 'Half the volume should give half the signal = %0.2f' % cut.getSignalArray()
+   print('Total signal in input = {}'.format(sum(signal)))
+   print('Half the volume should give half the signal = {}'.format(cut.getSignalArray()[0][0]))
 
 Output:
 
 .. testoutput:: ExampleMDHisto
 
-   Total signal in input = 100.00
-   Half the volume should give half the signal = 50.00
+   Total signal in input = 100.0
+   Half the volume should give half the signal = 50.0
    
 .. categories::
 
diff --git a/docs/source/algorithms/CylinderPaalmanPingsCorrection-v2.rst b/docs/source/algorithms/CylinderPaalmanPingsCorrection-v2.rst
index fa3ff1f0b80ded918b8f8d45d20d6037213e009e..fdc45aeb0b045aa0fcb88dc27766ebdb26a1a7c5 100644
--- a/docs/source/algorithms/CylinderPaalmanPingsCorrection-v2.rst
+++ b/docs/source/algorithms/CylinderPaalmanPingsCorrection-v2.rst
@@ -74,7 +74,7 @@ Usage
                                           Emode='Indirect',
                                           Efixed=1.845)
 
-    print 'Correction workspaces: %s' % (', '.join(corr.getNames()))
+    print('Correction workspaces: {}'.format((', '.join(corr.getNames()))))
 
 Output:
 
diff --git a/docs/source/algorithms/DNSComputeDetEffCorrCoefs-v1.rst b/docs/source/algorithms/DNSComputeDetEffCorrCoefs-v1.rst
index 3b7b0c7c738147187b9a11a7b03992325debed7f..531751a4edb624bd15171d245527f48db72f4ad0 100644
--- a/docs/source/algorithms/DNSComputeDetEffCorrCoefs-v1.rst
+++ b/docs/source/algorithms/DNSComputeDetEffCorrCoefs-v1.rst
@@ -81,23 +81,25 @@ Usage
    # Calculate correction coefficients
    coefs = DNSComputeDetEffCorrCoefs([vana_sf, vana_nsf], [bkgr_sf, bkgr_nsf])
 
-   print "First 3 correction coefficients: "
+   print("First 3 correction coefficients: ")
    for i in range(3):
-        print round(coefs.readY(i),2)
+        print(round(coefs.readY(i),2))
 
-   print "Is first detector masked?", coefs.getInstrument().getDetector(1).isMasked()
+   print("Is first detector masked? {}".format(coefs.getInstrument().getDetector(1).isMasked()))
 
    # load sample data 
    rawdata = LoadDNSLegacy('oi196012pbi.d_dat', Normalization='duration', CoilCurrentsTable=curtable)
 
    # apply correction
    corrected_data = rawdata/coefs
-   print "First 3 corrected data points"
+   print("First 3 corrected data points")
    for i in range(3):
-        print round(corrected_data.readY(i),2)
+        print(round(corrected_data.readY(i),2))
 
 Output:
 
+.. code-block:: python
+
    First 3 correction coefficients:
 
    0.0
diff --git a/docs/source/algorithms/DNSFlippingRatioCorr-v1.rst b/docs/source/algorithms/DNSFlippingRatioCorr-v1.rst
index 6e02d9065c02b8d0619b1870e9b2da3ca720c257..438f785749e0f19ba84182699065993067f9a08f 100644
--- a/docs/source/algorithms/DNSFlippingRatioCorr-v1.rst
+++ b/docs/source/algorithms/DNSFlippingRatioCorr-v1.rst
@@ -110,7 +110,7 @@ Usage
     vana_ratio = sf_corrected/nsf_corrected
 
     # ratio must be around 2, print first 5 points of the data array
-    print np.around(vana_ratio.extractY()[:5])
+    print(np.around(vana_ratio.extractY()[:5]))
 
 Output:
 
diff --git a/docs/source/algorithms/DNSMergeRuns-v1.rst b/docs/source/algorithms/DNSMergeRuns-v1.rst
index 1749f53b570ab7f6f7438b2044b35cf1f420c417..5221c1f2fea61f1e424f3a5868bec6ea6c1cf5b5 100644
--- a/docs/source/algorithms/DNSMergeRuns-v1.rst
+++ b/docs/source/algorithms/DNSMergeRuns-v1.rst
@@ -81,10 +81,10 @@ Usage
     for f in datafiles:
         try:
             wname = splitext(f)[0]
-            #print "Processing ", wname  # uncomment if needed
+            #print("Processing {}".format(wname))  # uncomment if needed
             LoadDNSLegacy(Filename=join(mypath, f), OutputWorkspace=wname, CoilCurrentsTable=coilcurrents, Normalization='duration')
         except RuntimeError as err:
-            print err
+            print(err)
         else:
             wslist.append(wname)
 
@@ -95,19 +95,19 @@ Usage
 
     # print selected values from merged workspaces
     two_theta = merged.extractX()[0]
-    print "First 5 2Theta values: ", two_theta[:5]
+    print("First 5 2Theta values: {}".format(two_theta[:5]))
     q = mergedQ.extractX()[0]
-    print "First 5 |Q| values: ", np.round(q[:5], 3)
+    print("First 5 |Q| values: {}".format(np.round(q[:5], 3)))
     d = mergedD.extractX()[0]
-    print "First 5 d values: ", np.round(d[:5], 3)
+    print("First 5 d values: {}".format(np.round(d[:5], 3)))
 
 Output:
 
-   First 5 2Theta values:  [ 7.5  8.   8.5  9.   9.5]
+   First 5 2Theta values: [ 7.5  8.   8.5  9.   9.5]
    
-   First 5 Q values:  [ 0.249  0.266  0.282  0.299  0.315]
+   First 5 Q values: [ 0.249  0.266  0.282  0.299  0.315]
    
-   First 5 d values:  [ 1.844  1.848  1.852  1.856  1.86 ]
+   First 5 d values: [ 1.844  1.848  1.852  1.856  1.86 ]
 
 .. categories::