diff --git a/Framework/PythonInterface/test/python/mantid/SimpleAPIFitTest.py b/Framework/PythonInterface/test/python/mantid/SimpleAPIFitTest.py
index 0bf4d83094136e4007de3507ec60f2615e952e01..a55a1a8852c5d0a31071356800a125996ab606b3 100644
--- a/Framework/PythonInterface/test/python/mantid/SimpleAPIFitTest.py
+++ b/Framework/PythonInterface/test/python/mantid/SimpleAPIFitTest.py
@@ -91,7 +91,6 @@ class SimpleAPIFitTest(unittest.TestCase):
         self.assertTrue(isinstance(retvals.OutputWorkspace, MatrixWorkspace))
         self.assertTrue(isinstance(retvals.FunctionString, str))
 
-
     def test_that_dialog_call_raises_runtime_error(self):
         try:
             FitDialog()
diff --git a/docs/source/algorithms/ExtractFFTSpectrum-v1.rst b/docs/source/algorithms/ExtractFFTSpectrum-v1.rst
index a1509a1f2091277bcb8b1abe72fdd2b4e8be8169..7d75cd5a0624cdeb8ea37030c4294d5cdfdf0c9e 100644
--- a/docs/source/algorithms/ExtractFFTSpectrum-v1.rst
+++ b/docs/source/algorithms/ExtractFFTSpectrum-v1.rst
@@ -57,6 +57,7 @@ Usage
 
 .. testcode:: Ex
 
+    from __future__ import print_function
     import numpy
 
     # Funtions x and y defined over the time domain: z(t) = x(t) + i * y(t)
@@ -79,12 +80,14 @@ Usage
 
     # Test the real part with a fitting to the expected Lorentzian
     myFunc = 'name=Lorentzian,Amplitude=0.05,PeakCentre=0,FWHM=6,ties=(PeakCentre=0)'
-    fitStatus, chiSq, covarianceTable, paramTable, fitWorkspace = Fit(Function=myFunc, InputWorkspace='wsfr', StartX=-40, EndX=40, CreateOutput=1)
+    fit_output = Fit(Function=myFunc, InputWorkspace='wsfr', StartX=-40, EndX=40, CreateOutput=1)
+    paramTable = fit_output.OutputParameters  # table containing the optimal fit parameters
     print("Theoretical FWHM = 1/(pi*tau)=6.367 -- Fitted FWHM value is: %.3f" % paramTable.column(1)[2])
 
     # Test the imaginary part with a fitting to the expected Gaussian
     myFunc = 'name=Gaussian,Height=0.1,PeakCentre=0,Sigma=3.0, ties=(PeakCentre=0)'
-    fitStatus, chiSq, covarianceTable, paramTable, fitWorkspace = Fit(Function=myFunc, InputWorkspace='wsfi', StartX=-15, EndX=15, CreateOutput=1)
+    fit_output = Fit(Function=myFunc, InputWorkspace='wsfi', StartX=-15, EndX=15, CreateOutput=1)
+    paramTable = fit_output.OutputParameters
     print("Theoretical Sigma = 1/(2*pi*sigma)=3.183 -- Fitted Sigma value is: %.3f" % paramTable.column(1)[2])
 
 Output:
diff --git a/docs/source/algorithms/Fit-v1.rst b/docs/source/algorithms/Fit-v1.rst
index 481b5c7d127f401daf42680ffbd2ef52abf0c582..e8cafc709b5b69bc70ee8d06b66280cecac3f4c3 100644
--- a/docs/source/algorithms/Fit-v1.rst
+++ b/docs/source/algorithms/Fit-v1.rst
@@ -382,6 +382,7 @@ Usage
 
 .. testcode:: ExFitPeak
 
+    from __future__ import print_function
    # create a workspace with a gaussian peak sitting on top of a linear (here flat) background
    ws = CreateSampleWorkspace(Function="User Defined", UserDefinedFunction="name=LinearBackground, \
       A0=0.3;name=Gaussian, PeakCentre=5, Height=10, Sigma=0.7", NumBanks=1, BankPixelWidth=1, XMin=0, XMax=10, BinWidth=0.1)
@@ -401,16 +402,17 @@ Usage
    #myFunc = 'name=LinearBackground, A0=0.3;name=Gaussian, Height='+height+', PeakCentre='+tryCentre+', Sigma='+sigma
 
    # Do the fitting
-   fitStatus, chiSq, covarianceTable, paramTable, fitWorkspace = Fit(InputWorkspace='ws', \
-      WorkspaceIndex=0, StartX = startX, EndX=endX, Output='fit', Function=myFunc)
+   fit_output = Fit(InputWorkspace='ws', WorkspaceIndex=0, StartX = startX, EndX=endX, Output='fit', Function=myFunc)
+   paramTable = fit_output.OutputParameters  # table containing the optimal fit parameters
+   fitWorkspace = fit_output.OutputWorkspace
 
-   print "The fit was: " + fitStatus
-   print("chi-squared of fit is: %.2f" % chiSq)
+   print( "The fit was: " + fit_output.OutputStatus)
+   print("chi-squared of fit is: %.2f" % fit_output.OutputChi2overDoF)
    print("Fitted Height value is: %.2f" % paramTable.column(1)[0])
    print("Fitted centre value is: %.2f" % paramTable.column(1)[1])
    print("Fitted sigma value is: %.2f" % paramTable.column(1)[2])
    # fitWorkspace contains the data, the calculated and the difference patterns
-   print "Number of spectra in fitWorkspace is: " +  str(fitWorkspace.getNumberHistograms())
+   print("Number of spectra in fitWorkspace is: " +  str(fitWorkspace.getNumberHistograms()))
    print("The 20th y-value of the calculated pattern: %.4f" % fitWorkspace.readY(1)[19])
 
 Output:
@@ -429,31 +431,32 @@ Output:
 
 .. testcode:: simFit
 
-	import math
-	import numpy as np
-
-	# create data
-	xData=np.linspace(start=0,stop=10,num=22)
-	yData=[]
-	for x in xData:
-		yData.append(2.0)
-	yData2=[]
-	for x in xData:
-		yData2.append(5.0)
-	# create workspaces
-	input = CreateWorkspace(xData,yData)
-	input2 = CreateWorkspace(xData,yData2)
-	# create function
-	myFunc=';name=FlatBackground,$domains=i,A0=0'
-	multiFunc='composite=MultiDomainFunction,NumDeriv=1'+myFunc+myFunc+";"
-	# do fit
-	fitStatus, chiSq, covarianceTable, paramTable, fitWorkspace = Fit( Function=multiFunc,\
-				InputWorkspace=input, WorkspaceIndex=0, \
-				InputWorkspace_1=input2, WorkspaceIndex_1=0, \
-				StartX = 0.1, EndX=9.5, StartX_1 = 0.1, EndX_1=9.5,Output='fit' )
-	# print results	
-	print "Constant 1: {0:.2f}".format(paramTable.column(1)[0])
-	print "Constant 2: {0:.2f}".format(paramTable.column(1)[1])
+    from __future__ import print_function
+    import math
+    import numpy as np
+
+    # create data
+    xData=np.linspace(start=0,stop=10,num=22)
+    yData=[]
+    for x in xData:
+        yData.append(2.0)
+    yData2=[]
+    for x in xData:
+        yData2.append(5.0)
+    # create workspaces
+    input = CreateWorkspace(xData,yData)
+    input2 = CreateWorkspace(xData,yData2)
+    # create function
+    myFunc=';name=FlatBackground,$domains=i,A0=0'
+    multiFunc='composite=MultiDomainFunction,NumDeriv=1'+myFunc+myFunc+";"
+    # do fit
+    fit_output = Fit(Function=multiFunc, InputWorkspace=input, WorkspaceIndex=0, \
+                     InputWorkspace_1=input2, WorkspaceIndex_1=0, \
+                     StartX = 0.1, EndX=9.5, StartX_1 = 0.1, EndX_1=9.5,Output='fit' )
+    paramTable = fit_output.OutputParameters  # table containing the optimal fit parameters
+    # print results
+    print("Constant 1: {0:.2f}".format(paramTable.column(1)[0]))
+    print("Constant 2: {0:.2f}".format(paramTable.column(1)[1]))
 
 
 Output:
diff --git a/docs/source/algorithms/LoadSassena-v1.rst b/docs/source/algorithms/LoadSassena-v1.rst
index c66d4566ff2e3d6175a7808936288e869acb7d5d..6209e30f95c515a5951bddb66ab72cc80e6ce8e8 100644
--- a/docs/source/algorithms/LoadSassena-v1.rst
+++ b/docs/source/algorithms/LoadSassena-v1.rst
@@ -54,6 +54,7 @@ Usage
 
 .. testcode:: Ex
 
+    from __future__ import print_function
     ws = LoadSassena("loadSassenaExample.h5", TimeUnit=1.0)
     print 'workspaces instantiated: ', ', '.join(ws.getNames())
     fqtReal = ws[1] # Real part of F(Q,t)
@@ -66,15 +67,16 @@ Usage
     myFunc = 'name=Gaussian,Height={0},PeakCentre={1},Sigma={2}'.format(intensity,center,sigma)
 
     # Call the Fit algorithm and perform the fit
-    fitStatus, chiSq, covarianceTable, paramTable, fitWorkspace =\
-    Fit(Function=myFunc, InputWorkspace=fqtReal, WorkspaceIndex=0, StartX = startX, EndX=endX, Output='fit')
+    fit_output = Fit(Function=myFunc, InputWorkspace=fqtReal, WorkspaceIndex=0, StartX = startX, EndX=endX, Output='fit')
+    paramTable = fit_output.OutputParameters  # table containing the optimal fit parameters
+    fitWorkspace = fit_output.OutputWorkspace
 
-    print "The fit was: " + fitStatus
+    print("The fit was: " + fit_output.OutputStatus)
     print("Fitted Height value is: %.2f" % paramTable.column(1)[0])
     print("Fitted centre value is: %.2f" % abs(paramTable.column(1)[1]))
     print("Fitted sigma value is: %.1f" % paramTable.column(1)[2])
     # fitWorkspace contains the data, the calculated and the difference patterns
-    print "Number of spectra in fitWorkspace is: " +  str(fitWorkspace.getNumberHistograms())
+    print( "Number of spectra in fitWorkspace is: " +  str(fitWorkspace.getNumberHistograms()))
     print("The 989th y-value of the fitted curve: %.3f" % fitWorkspace.readY(1)[989])
 
 Output:
diff --git a/docs/source/algorithms/SassenaFFT-v1.rst b/docs/source/algorithms/SassenaFFT-v1.rst
index 92b7f747f24229f85e9887117c98a09d90c490f7..fe370667e0b1f44f2aa7c0df96da2a16e7e993bf 100644
--- a/docs/source/algorithms/SassenaFFT-v1.rst
+++ b/docs/source/algorithms/SassenaFFT-v1.rst
@@ -66,6 +66,7 @@ Usage
 
 .. testcode:: Ex
 
+    from __future__ import print_function
     ws = LoadSassena("loadSassenaExample.h5", TimeUnit=1.0)
     SassenaFFT(ws, FFTonlyRealPart=1, Temp=1000, DetailedBalance=1)
 
@@ -82,15 +83,17 @@ Usage
     myFunc = 'name=Gaussian,Height={0},PeakCentre={1},Sigma={2}'.format(intensity,center,sigma)
 
     # Call the Fit algorithm and perform the fit
-    fitStatus, chiSq, covarianceTable, paramTable, fitWorkspace =\
-    Fit(Function=myFunc, InputWorkspace=sqt, WorkspaceIndex=0, StartX = startX, EndX=endX, Output='fit')
+    fit_output = Fit(Function=myFunc, InputWorkspace=sqt, WorkspaceIndex=0,
+                     StartX = startX, EndX=endX, Output='fit')
+    paramTable = fit_output.OutputParameters  # table containing the optimal fit parameters
+    fitWorkspace = fit_output.OutputWorkspace
 
-    print "The fit was: " + fitStatus
+    print("The fit was: " + fit_output.OutputStatus)
     print("Fitted Height value is: %.1f" % paramTable.column(1)[0])
     print("Fitted centre value is: %.1f" % abs(paramTable.column(1)[1]))
     print("Fitted sigma value is: %.4f" % paramTable.column(1)[2])
     # fitWorkspace contains the data, the calculated and the difference patterns
-    print "Number of spectra in fitWorkspace is: " +  str(fitWorkspace.getNumberHistograms())
+    print("Number of spectra in fitWorkspace is: " +  str(fitWorkspace.getNumberHistograms()))
 
 Output:
 
diff --git a/docs/source/concepts/Fitting.rst b/docs/source/concepts/Fitting.rst
index 45d8a6c237e17375785fda8e2f4c3fc180d3a689..1ecac7871de94cdc61ae485ed9a118f6d33c5393 100644
--- a/docs/source/concepts/Fitting.rst
+++ b/docs/source/concepts/Fitting.rst
@@ -159,6 +159,7 @@ implementation used for solving crystal structures from powder diffraction data.
     # For this demo example, just one fitting parameter is globally fitted, the peak center of a Gaussian peak
     # Please bear in mind the example here is to demonstrate this algorithm not provide a real global fitting problem
 
+    from __future__ import print_function
     from random import random
     from time import sleep
 
@@ -179,24 +180,25 @@ implementation used for solving crystal structures from powder diffraction data.
 
         # Do a fit from this starting value of the peak centre fitting parameter
         # Note choice of local minimizer will affect the outcome
-        d0, costFuncVal, d1, d2, d3 = Fit(InputWorkspace='data', WorkspaceIndex=0, \
-    	   StartX = startX, EndX=endX, Output='fit', \
-    	   Function='name=Gaussian,Height=10,PeakCentre='+tryCentre+',Sigma=20',
-    	   Minimizer='Conjugate gradient (Fletcher-Reeves imp.)')
+        fit_results = Fit(InputWorkspace='data', WorkspaceIndex=0, \
+                          StartX = startX, EndX=endX, Output='fit', \
+                          Function='name=Gaussian,Height=10,PeakCentre='+tryCentre+',Sigma=20',
+                          Minimizer='Conjugate gradient (Fletcher-Reeves imp.)')
+        costFuncVal = fit_results.OutputChi2overDoF
 
         # Here simply keep record of the best fit found, but this could easily be extended to
         # keep a record of all the minima found
         if costFuncVal < costFuncBest:
             costFuncBest = costFuncVal
             # here keep clone of best fit
-    	    CloneWorkspace(InputWorkspace='fit_Workspace',OutputWorkspace='fitBest')
+            CloneWorkspace(InputWorkspace='fit_Workspace',OutputWorkspace='fitBest')
 
         # Uncomment the sleep if would like to watch this algorithm trying to
         # find the global minima (graphically and/or from command line)
-        # print costFuncVal
+        # print(costFuncVal)
         # sleep(2)
 
-    print 'test'
+    print('test')
 
 .. testoutput:: LocalMinimizationRandowStartingPoints
     :hide:
diff --git a/docs/source/fitfunctions/TabulatedFunction.rst b/docs/source/fitfunctions/TabulatedFunction.rst
index 0b8239ab5d03d5b5e2517356d519c39c7a02f212..09771f6365b09300564df17c0f50e9c896b1a74c 100644
--- a/docs/source/fitfunctions/TabulatedFunction.rst
+++ b/docs/source/fitfunctions/TabulatedFunction.rst
@@ -32,6 +32,7 @@ Usage
 
 .. testcode:: Ex
 
+    from __future__ import print_function
     ws1=LoadNexus('tabulatedFunctionExample.nxs')
     
     # Clone the workspace by rescaling and shift 
@@ -41,14 +42,18 @@ Usage
     
     # Call the Fit algorithm and perform the fit
     myFunc='name=TabulatedFunction,Workspace=ws1,WorkspaceIndex=0,Scaling=1.0,Shift=0.0'
+    fit_output = Fit(Function=myFunc, InputWorkspace=ws2, Output='fit')
+
     fitStatus, chiSq, covarianceTable, paramTable, fitWorkspace =\
     Fit(Function=myFunc, InputWorkspace=ws2, Output='fit')
+    paramTable = fit_output.OutputParameters  # table containing the optimal fit parameters
+    fitWorkspace = fit_output.OutputWorkspace
 
-    print "The fit was: " + fitStatus
+    print("The fit was: " + fit_output.OutputStatus)
     print("Fitted Scaling value is: %.2f" % paramTable.column(1)[0])
     print("Fitted Shift value is: %.3f" % abs(paramTable.column(1)[1]))
     # fitWorkspace contains the data, the calculated and the difference patterns
-    print "Number of spectra in fitWorkspace is: " +  str(fitWorkspace.getNumberHistograms())
+    print("Number of spectra in fitWorkspace is: " +  str(fitWorkspace.getNumberHistograms()))
 
 Output:
 
diff --git a/docs/source/fitfunctions/TeixeiraWaterSQE.rst b/docs/source/fitfunctions/TeixeiraWaterSQE.rst
index 218a1e97ff976b8090583a76b1e9fbe5cc7b8445..df752d75a9ab2602fa358bc5cd57bfb87ed4b58b 100644
--- a/docs/source/fitfunctions/TeixeiraWaterSQE.rst
+++ b/docs/source/fitfunctions/TeixeiraWaterSQE.rst
@@ -68,16 +68,17 @@ The value of the momentum transfer :math:`Q` is contained in the loaded data
     name=LinearBackground"""
     # Let's fit spectrum with workspace index 5. Appropriate value of Q is picked up
     # automatically from workspace 'data' and passed on to the fit function
-    out,chi2,covariance,params,curves=Fit(Function=function, InputWorkspace=data,
-        WorkspaceIndex=5, CreateOutput=True, Output="fit", MaxIterations=100)
+    fit_output = Fit(Function=function, InputWorkspace=data, WorkspaceIndex=5,
+                     CreateOutput=True, Output="fit", MaxIterations=100)
+    params = fit_output.OutputParameters  # table containing the optimal fit parameters
+
     # Check some results
-    drow=params.row(6)
     DiffCoeff = params.row(6)["Value"]
     Tau = params.row(7)["Value"]
     if abs(DiffCoeff-2.1)/2.1 < 0.1 and abs(Tau-1.85)/1.85 < 0.1:
         print("Optimal parameters within 10% of expected values")
     else:
-        print(DiffCoeff, Tau, chi2)
+        print(DiffCoeff, Tau, fit_output.OutputChi2overDoF)
 
 
 **Example - Global fit to a synthetic signal:**