diff --git a/README.MD b/README.MD
index e05f64f..676ab7a 100644
--- a/README.MD
+++ b/README.MD
@@ -27,6 +27,10 @@ All Python packages that are needed in order to run the model are listed in the
 	
 	pip install -r requirements.txt
 	
+..then install the pathos framework and its modified externals from: https://github.com/uqfoundation
+A MPI distribution is needed for mpi4py which itself is a dependency of the pathos dependency pyina.
+
+The pathos framework is used for multiprocessing in steadyStateAnalysis.py only. Uncomment Line 73 to 88 and comment line 90 to 124 out to use a single core approach which does not require pathos.	
 It was noticed that under some systems LSODA integrator from the scipy.odeint package is not available. This integrator is necessary for a correct computation. Please check if it is available on your system before running the model.
 
 ## Example
diff --git a/requirements.txt b/requirements.txt
index dbc47ad..31a9441 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,10 +1,21 @@
-matplotlib==1.4.2
+dill==0.2.2
+klepto==0.1.1
+matplotlib==1.4.3
 mock==1.0.1
-nose==1.3.4
-numpy==1.9.1
+mpi4py==1.3.1
+mystic==0.2a2.dev0
+nose==1.3.6
+numpy==1.9.2
+pathos==0.2a1.dev0
+pox==0.2.1
+pp==1.6.4.4
+ppft==1.6.4.5
+processing===0.52-pathos
+pyina==0.2a1.dev0
 pyparsing==2.0.3
-python-dateutil==2.4.0
-pytz==2014.10
-scipy==0.15.0
+pyre===0.8.2.0-pathos
+python-dateutil==2.4.2
+pytz==2015.2
+scipy==0.15.1
 six==1.9.0
-wsgiref==0.1.2
+sympy==0.7.6
diff --git a/steadyStateAnalysis.py b/steadyStateAnalysis.py
index 63303a9..bae706a 100644
--- a/steadyStateAnalysis.py
+++ b/steadyStateAnalysis.py
@@ -5,7 +5,7 @@
 Calculate the steady state of the system when state transition are switched off
 
 
-Copyright (C) 2014-2015  Anna Matuszyńska, Oliver Ebenhöh
+Copyright (C) 2014-2015  Anna Matuszyńska, Oliver Ebenhöh, Philipp Norf
 
 This program is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
@@ -22,16 +22,16 @@
 """
 
 import numpy as np
-from scipy.integrate import ode
 
-import matplotlib.pyplot as plt
-from mpl_toolkits.mplot3d import Axes3D
+import sys
+import dill
+from functools import partial
 
 import petcModel
-import lightProtocol
 import parametersPETC
 import simulate
-from misc import pH, pHinv
+
+from time import time, sleep
 
 p = parametersPETC.ParametersPETC()
 # ---------------------------------------- #
@@ -44,61 +44,117 @@
 m = petcModel.PETCModel(p)
 s = simulate.Sim(m)
 
-PFDrange = np.linspace(75,1500,20)
-STrange = np.linspace(0,1,20)
+# Length of PFDrange and STrange set through a command line argument:
+if len(sys.argv) == 3 and sys.argv[2].isdigit() and int(sys.argv[2]) >= 3:
+    N = int(sys.argv[2])
+else:
+    N = 20 # default. This also prevent input errors.
 
-Ys = np.zeros([len(STrange), len(PFDrange), 8])
+PFDrange = np.linspace(75,1500,N)
+STrange = np.linspace(0,1,N)
 
-# dark adapted state. Not important
-y0=np.array([p.PQtot, 0.0202, 5.000, 0.0000, 0.0000, 0.0001, 0.9, 0.0000])
+Ys = np.zeros([N,N,8])
 
-for i in range(len(STrange)):
-	y0[6] = STrange[i]
-	print('Teraz liczymy dla ' + str(i)) 
-	Y = s.steadyStateLightScan(PFDrange,y0)
-	Ys[i,:] = Y
+# dark adapted state. Not important         #[[]] see line 51
+y0=np.array([[p.PQtot, 0.0202, 5.000, 0.0000, 0.0000, 0.0001, 0.9, 0.0000]])
 
-import pickle
-ss = {"STrange": STrange, "PFDrange": PFDrange, "Ys": Ys}
-output = open('steadyStateAnalysisFixedST.pkl', 'wb')
-pickle.dump(ss, output, 2)
-output.close()
+# This replaces the part of y0[6] = STrange[i] in the previous loop.
+y0 = y0.repeat(N,0) # this works because of the double brackets [[]] in y0
+y0[:,6] = STrange
 
+# Fix PFDrange parameter
+# This might be unnecassary, pathos.multiprocessing can deal with functions
+# which require multiple arguments. However, I haven't tested if it can deal
+# with different numbers of arguments, yet.
+steadyState = partial(s.steadyStateLightScan,PFDrange)
+
+# --------------------------------------- #
+#  This approach uses only a single core  #
+'''
+clock = time() # starting time
+for i in range(N):
+    print('Teraz liczymy dla ' + str(i)) # I don't get this line...
+    Y = steadyState(y0[i])
+    Ys[i,:] = Y
+clock = time()-clock # finishing time for N iterations
+
+# Serilisation of results and stats: 
+ss = {'STrange': STrange, 'PFDrange': PFDrange, 'Ys': Ys, 'Sec': clock}
+output = open('steadyStateAnalysisFixedST_SC_N' + str(N) + '.pkl', 'wb')
+dill.dump(ss,output,2)
+output.close()
 
+print(clock)
+'''
+# --------------------------------------- #
+#    This approach uses multiple cores    #
+from pathos.multiprocessing import ProcessingPool, cpu_count
+
+if __name__ == '__main__': # This is essential if Used on windows!
+    clock = time() # starting time
+    
+    # creates a worker pool from given comand line parameter. If the given
+    # parameter is to large all detectable CPUs will be utilised. If the given
+    # parameter is nonsense only 1 core will be utilized.
+    workers = 1
+    if len(sys.argv) >= 2 and sys.argv[1].isdigit() and int(sys.argv[1]) > 0:
+        workers = cpu_count()
+        if int(sys.argv[1]) <= workers:
+            workers = int(sys.argv[1])
+    
+    print 'N:  ' + str(N)
+    print 'PW: ' + str(workers)
+    sleep(3) # just 3 seconds pause to read the input again.
+
+    # All the magic happens here:
+    pool = ProcessingPool(workers)
+    Ys = pool.map(steadyState,y0)   
+
+    clock = time()-clock # elapsed time
+    print 'Seconds: ' + str(clock) # Not essential but useful.
+
+    # Serilisation of results and stats:
+    ss = {'STrange': STrange, 'PFDrange': PFDrange, 'Ys': Ys, 'Sec': clock, 'PoolWorkers': workers}
+    output = open('steadyStateAnalysisFixedST_MC_N' + str(N) + '.pkl', 'wb')
+    dill.dump(ss,output,2)
+    output.close()
+
+else:
+    print('Well, something went wrong.')
 
 #================================================================= #
 # 3 D plotting routine to obtain figure as in Ebenhoeh et al. 2014 #
 '''
-import pickle
+import dill
 from mpl_toolkits.mplot3d import Axes3D
-from matplotlib import cm
+import matplotlib
 from matplotlib.ticker import LinearLocator, FormatStrFormatter
 import matplotlib.pyplot as plt
 import numpy as np
 
-file = open('steadyStateAnalysisFixedST_corrected.pkl', 'rb')
-data = pickle.load(file)
+file = open('steadyStateAnalysisFixedST_MC_N20.pkl', 'rb')
+data = dill.load(file)
 
 
 ST = data['STrange']
 PFD = data['PFDrange']
 X,Y = np.meshgrid(PFD, ST)
 
-Yss = np.zeros([len(ST),20])
+Yss = np.zeros([len(ST),len(data['Ys'])])
 
 for i in range(len(ST)):
         for j in range(len(PFD)):
          Yss[i,j] = 1 - data['Ys'][i][j][0] / 17.5
 
-y_formatter = matplotlib.ticker.ScalarFormatter(useOffset = False)
-ax.yaxis.set_major_formatter(y_formatter)
-
 cm = matplotlib.cm.get_cmap('jet')
 fig = plt.figure()
 ax = fig.add_subplot(111, projection='3d')
 surf = ax.plot_surface(X,Y,Yss, rstride=1, cstride=1, cmap=cm,
                        linewidth=1, antialiased=True)
 
+y_formatter = matplotlib.ticker.ScalarFormatter(useOffset = False)
+ax.yaxis.set_major_formatter(y_formatter)
+
 ax.zaxis.set_major_locator(LinearLocator(10))
 ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
 
diff --git a/steadyStateAnalysisFixedST_MC_N100.pkl b/steadyStateAnalysisFixedST_MC_N100.pkl
new file mode 100644
index 0000000..5fbc71c
Binary files /dev/null and b/steadyStateAnalysisFixedST_MC_N100.pkl differ
diff --git a/steadyStateAnalysisFixedST_MC_N50.pkl b/steadyStateAnalysisFixedST_MC_N50.pkl
new file mode 100644
index 0000000..f71af55
Binary files /dev/null and b/steadyStateAnalysisFixedST_MC_N50.pkl differ