diff --git a/trunk/BD_MODELS/sOTConst.txt b/trunk/BD_MODELS/sOTConst.txt
new file mode 100644
index 0000000000000000000000000000000000000000..a98da38d3263464c3c3b0a838c1588159138c15f
--- /dev/null
+++ b/trunk/BD_MODELS/sOTConst.txt
@@ -0,0 +1,2 @@
+fax=1*p1
+fay=p2
diff --git a/trunk/BD_MODELS/sOTConst2.txt b/trunk/BD_MODELS/sOTConst2.txt
new file mode 100644
index 0000000000000000000000000000000000000000..42071d10f6545365b3637689a3416bbf6b760f27
--- /dev/null
+++ b/trunk/BD_MODELS/sOTConst2.txt
@@ -0,0 +1,2 @@
+fax=p3
+fay=p3
diff --git a/trunk/src/mse/DISPERSIONS/diffusionModel.py b/trunk/src/mse/DISPERSIONS/diffusionModel.py
index 8184738722cf35efeca7b36b38e8bfa46ee61d77..ef13327c5885de651bd9797d4e898b2b6b017d7d 100644
--- a/trunk/src/mse/DISPERSIONS/diffusionModel.py
+++ b/trunk/src/mse/DISPERSIONS/diffusionModel.py
@@ -1,15 +1,19 @@
-from mse.MODELS.models import model
+from mse.MODELS.models import model, strForDolf, formatCovariables
 import mse.dynamicalSystem
-from dolfinx.fem import Constant,form
+from dolfinx.fem import Constant,form,assemble_scalar
 from petsc4py.PETSc import ScalarType
 from dolfinx.fem.petsc import assemble_matrix
-from ufl import grad,dot
+from ufl import grad,dot, inner
 from mse.linAlgTools import matrix,vector
+from mse.toolsEnv import computerEnv, MSEPrint, MSEPrint_0
+from os.path import exists
+from sympy import diff
+import numpy as np
 
 class diffusionModel(model):
     def __init__(self,aDs,aDiffCoef=1) -> None:
         super().__init__(aDs)
-        self.coefName=aDiffCoef
+        aDs.addSecondOrderTerm(self)
         self.diffCoef=Constant(self.dS.dolfinMesh, ScalarType(aDiffCoef))
         self.block=None
         self.blockAdjoint=None
@@ -23,6 +27,8 @@ class diffusionModel(model):
     def _computeBlock(self):
         aux=dot(grad(self.diffCoef*self.dS.u),grad(self.dS.v))*self.dS.dx
         self.bilinearForm=form(aux)
+        print("aux = ",aux)
+        print("bilineaire form in diff model(cst): ",self.bilinearForm)
         self.block=matrix(assemble_matrix(self.bilinearForm))
     def updateBlockMat(self,compi,compj):
         if (compi==compj):
@@ -34,3 +40,171 @@ class diffusionModel(model):
     def updateBlockMatAdjoint(self,compi,compj):
         if (compi==compj):
             return self.blockAdjoint
+
+class txtDiffusionModel(model):
+    """ Class used to create a second order term with a txt file
+
+    """
+    unknownName="abcdefgh"
+    def _subFXForm(self,strModel,X):
+        for j in range(self.nbComps):
+            strModel=strModel.replace("$"+txtDiffusionModel.unknownName[j],"self.dS."+X+".comps["+str(j)+"]")
+        code="form(("+strModel+")"
+        model.myPrint("X ->"+code)
+        return code
+
+    def __init__(self,aDs, aTxtFile) -> None:
+        super().__init__(aDs)
+        aDs.addSecondOrderTerm(self)
+        self.txtFile=aTxtFile
+        self.block=None
+        self.blockAdjoint=None
+        self.bilinearForm=None
+        self._computeBlock()
+        self.updateLinearisedMat()
+    def setDiffCoef(self,aTxtFile):
+        self.diffCoef=Constant(self.dS.dolfinMesh, ScalarType(aDiffCoef))
+        self._computeBlock()
+        self.updateLinearisedMat()
+    def _computeBlock(self):
+        # we check if file exist
+        if (exists(self.txtFile+".txt")):
+            self.txtFile=self.txtFile
+        else:
+            self.txtFile=computerEnv.modelsDir+self.txtFile #self.txtFile form is "affineGrowth"
+        model.myPrint("Looking for model  "+self.txtFile)
+        try:
+            f = open(self.txtFile+".txt", "r")
+        except IOError:
+            model.myPrint("IOError in txtDiffModel when reading "+self.txtFile+".txt")
+            return
+        # we check if the derivate of parameter file exist
+        try:
+            fd = open(self.txtFile+"DVP.txt", "r")
+        except FileNotFoundError:
+            # if didn't exist, we create it
+            self.dvFileExist=False
+            print("Création fichier "+self.txtFile+"DVP.tx")
+            fd = open(self.txtFile+"DVP.txt", "w") #open file in mod write
+            # loop for each comp
+            for numComps in range(self.nbComps):
+                #loop for each dimension(=2)
+                for strDim in ['x','y']:
+                    strF=f.readline()
+                    strF=strF.split('=')[1].strip()#.replace("^","**")
+                    #loop for each param
+                    for p in range(self.dS.study.nbParam):
+                        grad1s=diff(strF,'p'+str(p+1))
+                        grad1=str(grad1s).replace("**","^")
+                        #write in fic
+                        fd.write("dp"+strDim+str(p+1)+"f"+txtDiffusionModel.unknownName[numComps]+"="+grad1+"\n")
+            f.close()
+            fd.close()
+            f = open(self.txtFile+".txt", "r")
+            try:
+                fd = open(self.txtFile+"DVP.txt", "r")
+            except IOError:
+                model.myPrint("IOError in txtDiffusionModel when reading "+self.txtFile+"DVP.txt")
+                return
+        except IOError:
+            model.myPrint("IOError in txtDiffusionModel when reading "+self.txtFile+"DVP.txt")
+            return
+        ### find idea for derivate in obprocess.computeGradV
+        ### TODO: read file + create varfs like txtmodel
+        nbDim = 2
+        self.varfs=[None]*(self.nbComps*(nbDim*self.dS.study.nbParam+1))#FIXME
+        #self.varfsTm1=[None]*self.nbComps
+        #for the comp number numComp \in {0,..,self.nbComps-1}:
+        #self.varfs[numComp*(self.nbComps+1)] contains the varf F_numComp for the rhs
+        #self.varfs[numComp*(self.nbComps+1)+(j+1)] contains the derivative d_jF_numComp, j \in {0,..,self.nbComps-1}
+        strkeys=r"[\s\+\-\*/\=\(\)\^]+"
+        for numComps in range(self.nbComps):
+            ##READ SECOND ORDER TERM
+            # Dxx
+            strlineX=f.readline().split('=')[1].strip()
+            strlineX=strForDolf(strlineX, self.unknownName)
+            strlineX=strlineX.replace('^','**')
+            # Dyy
+            strlineY=f.readline().split('=')[1].strip()
+            strlineY=strForDolf(strlineY, self.unknownName)
+            strlineY=strlineY.replace('^','**')
+            ##replace param
+            for count_p in range(len(self.dS.paramD)):
+                strlineX=strlineX.replace('p'+str(count_p+1),'self.dS.paramD['+str(count_p)+']')
+                strlineY=strlineY.replace('p'+str(count_p+1),'self.dS.paramD['+str(count_p)+']')
+            #### replace COVARIABLES
+            strlineX=formatCovariables(strlineX,"self.dS.dicCov[\"", "\"]")
+            strlineY=formatCovariables(strlineY,"self.dS.dicCov[\"", "\"]")
+            #### BUILD VARF USING FK
+            # v1: with Grad[i]
+            strline = strlineX + "*grad(self.dS.u)[0]*grad(self.dS.v)[0]*self.dS.dx+"+strlineY+"*grad(self.dS.u)[1]*grad(self.dS.v)[1]*self.dS.dx"
+            model.myPrint(strline)
+            print("varfs = ",strline)
+            code="self.varfs["+str(numComps)+"*(2*self.dS.study.nbParam+1)]="+self._subFXForm(strline,"uk")+")"
+            print("num index in varf =",numComps*(2*self.dS.study.nbParam+1))
+            print("code = ",code)
+            exec(compile(code, 'sumstring', 'exec'))
+                ##### BUILD VARF USING UTM1
+                #code="self.varfsTm1["+str(numComps)+"]="+self._subFXForm(strline,"utm1")+"*self.dS.v*self.dS.dx)"
+                #model.myPrint(code)
+                #exec(compile(code, 'sumstring', 'exec'))
+            # v2: with D= const(np.array[[Dxx,0],[0,Dyy]])
+            #Dxy=[None]*2
+            #code="Dxy[0]="+strlineX
+            #print("code Dxx = ",code)
+            #exec(compile(code, 'sumstring', 'exec'))
+            #code="Dxy[1]="+strlineY
+            #print("code Dyy = ",code)
+            #exec(compile(code, 'sumstring', 'exec'))
+            #print("Dxx = "+str(Dxy[0])+", Dyy = ",Dxy[1])
+            ### Mat
+            ##D = Constant(self.dS.dolfinMesh, np.array([[Dxy[0], 0], [0, Dxy[1]]], dtype=np.float64))
+            ### Vec
+            #D = Constant(self.dS.dolfinMesh, np.array([Dxy[0], Dxy[1]], dtype=np.float64))
+            ### R
+            ##D = Constant(self.dS.dolfinMesh, self.dS.study.param[0])
+            #strline = "dot(D*grad(self.dS.u)+ self.dS.u*grad(D),grad(self.dS.v))*self.dS.dx"
+            #code="self.varfs["+str(numComps)+"*(self.nbComps+1)]="+self._subFXForm(strline,"uk")+")"
+            #print("num index in varf =",numComps*(self.nbComps+1))
+            #print("code = ",code)
+            ##exec(compile(code, 'sumstring', 'exec'))
+            #self.varfs[0*(self.nbComps+1)]=form((dot(dot(D,grad(self.dS.u))+ dot(self.dS.u,grad(D)),grad(self.dS.v))*self.dS.dx))
+            #
+            ## READ DERIVATIVES PARAM
+            # NOTE: the value didn't depend of 'u', so didn't need 'self.subxForm'
+            # NOTE: we will use 'lambda' to create function ref to param.
+            # loop for each dim (=2)
+            for numDim in range(nbDim):
+                # loop for each param
+                for countP in range(len(self.dS.paramD)):
+                    strline=fd.readline().split('=')[1].strip()
+                    print("num index for derivate param =",numComps*(2*self.dS.study.nbParam+1)+countP+1+(numDim)*(self.dS.study.nbParam))
+                    print("strline==\"0\":",strline=="0")
+                    if strline == '0':
+                        pass
+                    else:
+                        strline=strForDolf(strline,self.unknownName)
+                        ##replace param
+                        for tmpCountP,tmpP in enumerate(self.dS.paramD):
+                            strline=strline.replace('p'+str(tmpCountP+1),'self.dS.study.param['+str(tmpCountP)+']')
+                        #replace COVARIABLES
+                            #tmpDPSOT = assemble_scalar(form((tmpDPSOTX + tmpDPSOTY)))
+                        strline=formatCovariables(strline,"self.dS.dicCov[\"", "\"]")
+                        model.myPrint("txtDiffusionModel: read txt model d"+txtDiffusionModel.unknownName[numComps]+"F"+txtDiffusionModel.unknownName[numComps]+"="+strline)
+                        code="self.varfs["+str(numComps)+"*(2*self.dS.study.nbParam+1)+"+str(tmpCountP)+"+1+"+str(numDim)+"*(self.dS.study.nbParam)]=lambda: "+strline
+                        print("codeDvp = ",code)
+                        exec(compile(code, 'sumstring', 'exec'))
+            for i in self.varfs:print (i)
+    def updateBlockMat(self,compi,compj):
+        if (compi==compj):
+            #return matrix(assemble_matrix(self.varfs[compi*(self.nbComps+1)]))
+            return matrix(assemble_matrix(self.varfs[0]))
+    def updateBlockMatAdjoint(self,compi,compj):
+        if (compi==compj):
+            return self.blockAdjoint
+    ### surcharge function "getlineraized" to add update
+    def getLinearisedMat(self):
+        super().updateLinearisedMat()
+        return self.mat
+        #NOTE: super()getLinearizedMat didn't work....
+    # def get(): super().update, super.get()
diff --git a/trunk/src/mse/obsProcess.py b/trunk/src/mse/obsProcess.py
index 5b3e8dd1cfb43077e5709130354c484e9c972343..622e60506efd116ee39d403244541763ef3a4c61 100644
--- a/trunk/src/mse/obsProcess.py
+++ b/trunk/src/mse/obsProcess.py
@@ -14,6 +14,7 @@ from mse.systemAdjoint import systemAdjoint
 from mse.TIMESCHEMES.timeScheme import implicitLinearBackwardTimeScheme
 from mse.timeDoOneStep import timeDoOneStepBackward
 from mse.simulator import simulatorBackward
+from mse.toolsEnv import MSEPrint, MSEPrint_0  
 from os.path import exists
 from os import mkdir
 class obsProcess:
@@ -23,6 +24,8 @@ class obsProcess:
     """
     def __init__(self, aDS, aDirURef=None)->None:
         #Dynamical system
+        #self.myPrint=MSEPrint_0
+        self.myPrint=MSEPrint
         self.dS=aDS
         self.dS.setObsProcess(self)
         self.index = self.dS.indexInSystem
@@ -78,7 +81,7 @@ class obsProcess:
             try:
                 f = open(sourceTerm.txtFile+".txt", "r")
             except IOError:
-                model.myPrint("IOError in txtModel when reading "+self.txtFile+".txt")
+                self.myPrint("IOError in txtModel when reading "+self.txtFile+".txt")
                 return 
             try:
                 fp = open(sourceTerm.txtFile+"Dp"+".txt", "r")
@@ -96,7 +99,7 @@ class obsProcess:
                         #write in fic
                         fp.write("d"+"p"+str(k+1)+"f"+sourceTerm.unknownName[numComp]+"="+grad1+"\n")
             except IOError:
-                model.myPrint("IOError in obsProcess when reading "+sourceTerm.txtFile+"Dp.txt")
+                self.myPrint("IOError in obsProcess when reading "+sourceTerm.txtFile+"Dp.txt")
                 return 
             f.close()
             fp.close()
@@ -117,7 +120,7 @@ class obsProcess:
 #            try:
 #                fp = open(sourceTerm.txtFile+"Dp"+".txt", "r")
 #            except IOError:
-#                model.myPrint("IOError in obsProcess when reading "+sourceTerm.txtFile+"Dp.txt")
+#                self.myPrint("IOError in obsProcess when reading "+sourceTerm.txtFile+"Dp.txt")
 #                return 
 #            for numComp in range(self.dS.nbComps):
 #               for numParam in range(self.dS.study.nbParam):
@@ -148,7 +151,7 @@ class obsProcess:
             try:
                 fp = open(sourceTerm.txtFile+"Dp"+".txt", "r")
             except IOError:
-                model.myPrint("IOError in obsProcess when reading "+sourceTerm.txtFile+"Dp.txt")
+                self.myPrint("IOError in obsProcess when reading "+sourceTerm.txtFile+"Dp.txt")
                 return 
             for numComp in range(self.dS.nbComps):
                 for numParam in range(self.dS.study.nbParam):
@@ -173,7 +176,10 @@ class obsProcess:
         pastT=self.dS.study.simulatorBackward.t0
         pastF=np.zeros(self.dS.study.nbParam)
         curF =np.zeros(self.dS.study.nbParam)
+        ### todo: READ SECONDoRDERtER.varfs[??], if != None add in grad
         # TODO: reprendre comme model: tableau de ufl, 1 dim de taille nbparam*???
+        tmpDPSOTX = 0
+        tmpDPSOTY = 0
         while(l_aExplorTrajAdjoint.replay()):
             #adjoint is load in ds.utm1 after call of l_aExplorTrajAdjoint.replay
             # coef si integration rectangle ou trapeze
@@ -199,7 +205,35 @@ class obsProcess:
                         #read u_a
                         code="self.daF.x.array[:]="+'+'.join(strdFp[i])
                         exec(compile(code,'sumstring', 'exec'))
-                        curF[i] = assemble_scalar(form(comp*self.daF*self.dS.dx))
+                        # 2 DIM
+                        # derivate of 2ndOrderTerm
+                        # X
+                        if (not hasattr(self.dS.secondOrderTerm,"varfs")):tmpDPSOT=0
+                        else:
+                            indX = count*(2*self.dS.study.nbParam+1)+i+1
+                            indY = count*(2*self.dS.study.nbParam+1)+i+1+self.dS.study.nbParam
+                            DpDx = self.dS.secondOrderTerm.varfs[indX]
+                            DpDy = self.dS.secondOrderTerm.varfs[indY]
+                            if(DpDx == None):tmpDPSOTX = 0*self.dS.dx
+                            else:
+                                tmpDPSOTX = DpDx()*grad(self.dS.trajState.comps[count])[0]*grad(comp)[0]*self.dS.dx
+                            # Y
+                            if(DpDy == None):tmpDPSOTY = 0*self.dS.dx
+                            else:
+                                tmpDPSOTY = DpDy()*grad(self.dS.trajState.comps[count])[1]*grad(comp)[1]*self.dS.dx
+                            #print("type of tmpDPSOTX: ",type(tmpDPSOTX))
+                            #print("(tmpDPSOTX): ",tmpDPSOTX)
+                            #print("type of form(tmpDPSOTX): ",type(form(tmpDPSOTX)))
+                            #print("type of :tmpDPSOTY ",type(tmpDPSOTY))
+                            #print("(tmpDPSOTY): ",tmpDPSOTY)
+                            #print("type of form(tmpDPSOTY): ",type(form(tmpDPSOTY)))
+                            #print("type of :form(comp*self.daF*self.dS.dx) ",type(form(comp*self.daF*self.dS.dx)))
+                            #print("type of : form((tmpDPSOTX + tmpDPSOTY))) ",type(form((tmpDPSOTX + tmpDPSOTY))))
+                            #print("type of : assemble_scalar(form((tmpDPSOTX + tmpDPSOTY)))) ",type(assemble_scalar(form((tmpDPSOTX + tmpDPSOTY)))))
+                            #print("assemble_scalar(form((tmpDPSOTX + tmpDPSOTY)))) ",assemble_scalar(form((tmpDPSOTX + tmpDPSOTY))))
+                            tmpDPSOT = assemble_scalar(form(tmpDPSOTX + tmpDPSOTY))
+                            print("tmpDPSOT: ",tmpDPSOT)
+                        curF[i] = assemble_scalar(form(comp*self.daF*self.dS.dx)) + tmpDPSOT
                         ### End V1
                         ### V2
                         # NOTE: cette version demande de creer un vecteur. L'impementation est à l'abandon
diff --git a/trunk/tests/test_StudyParam2.py b/trunk/tests/test_StudyParam2.py
index 5bedb024d3da5e5f92628116e045683d680f65a1..944658aa47a8fe3cd04e7975395589cf475ecd40 100644
--- a/trunk/tests/test_StudyParam2.py
+++ b/trunk/tests/test_StudyParam2.py
@@ -80,7 +80,3 @@ def test_Param2(result, index, param, trueParam):
     print("final mass = ",mass[0])
     assert r
 
-
-
-
-
diff --git a/trunk/tests/test_backwardModel.py b/trunk/tests/test_backwardModel.py
index 9e6af3619f847a54faeafd70c8d4795e78c5c7f0..511c2df3a41bd3a976228f78a306eb34be813fb7 100644
--- a/trunk/tests/test_backwardModel.py
+++ b/trunk/tests/test_backwardModel.py
@@ -90,7 +90,7 @@ quaDiff=quadraticDiff(ds2d1)
 
 # We test with different epsilon and parameters.
 print("Param a_ref =",aStudy.getParam)
-@pytest.mark.parametrize("defaultParam",[([0.5,0.1]),([1.5,1.])])
+@pytest.mark.parametrize("defaultParam",[([0.5,0.1]),([1.5,1.]),([1.2,.6])])
 def test_GradWithFiniteDifference(defaultParam):
     """
     test the approximation of GradV mse with finite Difference.
diff --git a/trunk/tests/test_secndOrdrTermTxt.py b/trunk/tests/test_secndOrdrTermTxt.py
new file mode 100644
index 0000000000000000000000000000000000000000..b4cb2b81783a80c78cb861d5ff470d5bdd8dd813
--- /dev/null
+++ b/trunk/tests/test_secndOrdrTermTxt.py
@@ -0,0 +1,90 @@
+from mse.study import study
+from mse.dynamicalSystem import dynamicalSystem,computeMass
+from mse.DISPERSIONS.diffusionModel import diffusionModel, txtDiffusionModel
+from mse.MODELS.models import txtModel
+from mse.system import system
+from mse.simulator import simulator
+from mse.timeDoOneStep import timeDoOneStep
+from mse.TIMESCHEMES.timeScheme import implicitLinearTimeScheme,implicitNonLinearTimeScheme
+from mse.toolsEnv import computerEnv,MSEPrint
+from mse.explorTrajectory import explorTrajectory
+
+import numpy as np
+import pytest
+from os import  mkdir
+from os.path import isdir
+
+
+
+### NOTE: without sourceTerm ? only 2nd order term read in file ?
+## test with constant sourceTerm + other.
+
+
+#un modele
+nameModele="affineGrowth"
+nameSendOrdTerm="sOTConst"
+sourceModele = computerEnv.modelsDir+nameModele+".txt"
+nbComps=1
+WITH_AFFINE_GROWTH=1
+nbparam=2
+eps=1e-6
+eps2=1e-2
+testNLTS=0
+
+    #une etude
+studName="TEST_SNDORDTERM"
+aStudy=study(studName,nbparam)
+aStudy.setGeometry("SQUARE",0.4) 
+    #un systeme
+aSystem=system(aStudy)
+    #un systeme dynamique
+ds2d1=dynamicalSystem(aStudy,2,"surface1.xdmf",nbComps=nbComps)
+     
+aStudy.setParam([0.,1.])
+aTS=None
+#schema en temps implicite lineaire
+aTS=implicitLinearTimeScheme()
+aTS.computeMatOnce=1
+    
+simulator(0,1,aStudy)
+aTStep=timeDoOneStep(aStudy,aTS,0.1)
+mass=np.zeros((ds2d1.nbComps),dtype=np.double)
+
+#add diffusion motion
+aDiffModel=txtDiffusionModel(ds2d1, nameSendOrdTerm)
+#aDiffModel=diffusionModel(ds2d1)
+#aDiffModel.setDiffCoef(1.)
+if (WITH_AFFINE_GROWTH):
+    m=txtModel(ds2d1,nameModele)
+    m.computeMatOnce=1
+    m.computeRhsOnce=1
+
+aTS=None
+# select the time algorithm
+if (testNLTS):
+    aTS=implicitNonLinearTimeScheme()
+else:
+    aTS=implicitLinearTimeScheme()
+#simu in t \in [0,1]
+simulator(0,1,aStudy)
+#based on fix time step 0.05
+aTStep=timeDoOneStep(aStudy,aTS,0.05)
+#run the simulation
+aStudy.simulator.doSimu()
+
+print("param = ",aStudy.param)
+#for visu, export simulation steps to vtk
+aStudy.simulator.exportSimuToVtk()
+
+#loop on simu to compute mass
+import numpy as np
+mass=np.zeros((ds2d1.nbComps),dtype=np.double)
+aExplorTraj=explorTrajectory(aStudy,aStudy.simuDir)
+aExplorTraj.rewind()
+while(aExplorTraj.replay()):
+    computeMass(ds2d1,mass)
+    print("mass "+str(aExplorTraj.curStep)+" = ",end="")
+    for m in mass:
+        print(m,end=", ")
+    print("")
+aStudy.end()