diff --git a/Compiler/BackEnd/BackendDAEUtil.mo b/Compiler/BackEnd/BackendDAEUtil.mo
index 7d5defd1b7e..894de16c1ff 100644
--- a/Compiler/BackEnd/BackendDAEUtil.mo
+++ b/Compiler/BackEnd/BackendDAEUtil.mo
@@ -3764,6 +3764,21 @@ algorithm
 end match;
 end traverseStmtsExpList;
 
+public function incidenceMatrixToSparsePattern
+" Applies absolute value to all entries in the incidence matrix and
+  also sort and unique them."
+  input BackendDAE.IncidenceMatrix m;
+  output BackendDAE.IncidenceMatrix res;
+protected
+  list<list<Integer>> lst,lst_1;
+algorithm
+  lst := arrayList(m);
+  lst_1 := List.mapList(lst, intAbs);
+  lst_1 := List.map1(lst_1, List.sort, intGt);
+  lst_1 := List.map1(lst_1, List.sortedUnique, intGe);
+  res := listArray(lst_1);
+end incidenceMatrixToSparsePattern;
+
 /******************************************************************
  stuff to calculate enhanced Adjacency matrix
 
diff --git a/Compiler/BackEnd/SymbolicJacobian.mo b/Compiler/BackEnd/SymbolicJacobian.mo
index 4737f041eea..0d4cddbbb82 100644
--- a/Compiler/BackEnd/SymbolicJacobian.mo
+++ b/Compiler/BackEnd/SymbolicJacobian.mo
@@ -1485,7 +1485,7 @@ algorithm
   end matchcontinue;
 end getSparsePatternHelp2;
 
-protected function transposeSparsePattern
+public function transposeSparsePattern
   input list<list<Integer>> inSparsePattern;
   input array<list<Integer>> inAccumList;
   input Integer inValue;
@@ -1503,6 +1503,37 @@ algorithm
   end for;
 end transposeSparsePattern;
 
+public function transposeSparsePatternTuple
+  input list<tuple<Integer, list<Integer>>> inSparsePattern;
+  input array<tuple<Integer,list<Integer>>> inAccumList;
+  output array<tuple<Integer,list<Integer>>> outSparsePattern = inAccumList;
+protected
+  Integer value;
+  list<Integer> tmplist;
+  list<Integer> oneList;
+  tuple<Integer,list<Integer>> tmpTuple;
+  Integer i;
+algorithm
+  for oneListTuple in inSparsePattern loop
+    (value, oneList) := oneListTuple;
+    for oneElem in oneList loop
+      tmpTuple := arrayGet(outSparsePattern,oneElem+1);
+      (_, tmplist) := tmpTuple;
+      tmplist := value::tmplist;
+      tmpTuple := (oneElem, tmplist);
+      MetaModelica.Dangerous.arrayUpdateNoBoundsChecking(outSparsePattern, oneElem+1, tmpTuple);
+    end for;
+  end for;
+  // sort all transposed lists
+  for i in 1:listLength(inSparsePattern) loop
+    tmpTuple := arrayGet(outSparsePattern,i);
+    (value, tmplist) := tmpTuple;
+    tmplist := List.sort(tmplist, intGt);
+    tmpTuple := (value, tmplist);
+    MetaModelica.Dangerous.arrayUpdateNoBoundsChecking(outSparsePattern, i, tmpTuple);
+  end for;
+end transposeSparsePatternTuple;
+
 protected function mapIndexColors
   input array<Integer> inColors;
   input Integer inMaxIndex;
diff --git a/Compiler/SimCode/SimCodeUtil.mo b/Compiler/SimCode/SimCodeUtil.mo
index 880a0c3a631..14ad7c53ff6 100644
--- a/Compiler/SimCode/SimCodeUtil.mo
+++ b/Compiler/SimCode/SimCodeUtil.mo
@@ -438,7 +438,7 @@ algorithm
     if Flags.getConfigBool(Flags.DAE_MODE) then
       daeModeSP := createDaeModeSparsePattern(bdaeModeVars, bdaeModeEqns, shared, crefToSimVarHT);
       residualVars := rewriteIndex(residualVars, 0);
-      algebraicVars := rewriteIndex(algebraicVars, 0);
+      algebraicVars := sortSimVarsAndWriteIndex(algebraicVars, crefToSimVarHT);
       daeModeData := SOME(SimCode.DAEMODEDATA(daeEquations, daeModeSP, residualVars, algebraicVars));
     else
       daeModeData := NONE();
@@ -1764,6 +1764,8 @@ algorithm
   foldArg := ({}, {}, {}, iuniqueEqIndex, itempvars, outDAEVars, outDAEEqns);
   (odaeEquations, oresidualVars, oalgebraicVars, ouniqueEqIndex, otempvars, outDAEVars, outDAEEqns) :=
     List.fold1(inSysts, createDAEEqsPrepare, shared, foldArg);
+  oresidualVars := listReverse(oresidualVars);
+  oalgebraicVars := listReverse(oalgebraicVars);
 end createDAEEquations;
 
 protected function createDAEEqsPrepare
@@ -1978,11 +1980,26 @@ protected
   BackendDAE.SparseColoring coloring;
   SimCode.JacobianMatrix simCodeJac;
   Integer nonZeroElements;
+  list<SimCodeVar.SimVar> inSimVars;
+  list<tuple<SimCodeVar.SimVar, BackendDAE.Var>> tupleVars;
+  list<BackendDAE.Var> inVarsLst;
+  BackendDAE.Variables sortedVars;
+  SimCodeVar.SimVar svar; BackendDAE.Var bvar;
 algorithm
   try
+    // sort variables in SimCode order
+    inVarsLst := BackendVariable.varList(inVars);
+    varCrefsList := listReverse(List.map(inVarsLst, BackendVariable.varCref));
+    inSimVars := getSimVars2Crefs(varCrefsList, crefToSimVarHT);
+    tupleVars := List.threadTuple(inSimVars, inVarsLst);
+    tupleVars := List.sort(tupleVars, compareSimVarTupleIndexGt);
+    inVarsLst := List.unzipSecond(tupleVars);
+    sortedVars := BackendVariable.listVar(listReverse(inVarsLst));
+    varCrefsList := List.map(inVarsLst, BackendVariable.varCref);
+
     // create adjacency matrix
     funcs := BackendDAEUtil.getFunctions(inShared);
-    system := BackendDAEUtil.createEqSystem(inVars, inEqns);
+    system := BackendDAEUtil.createEqSystem(sortedVars, inEqns);
     (system,_,_,_,_) := BackendDAEUtil.getIncidenceMatrixScalar(system, BackendDAE.SOLVABLE(), SOME(funcs));
     if debug then
       print("createDaeModeSparsePattern DAE-System:\n");
@@ -1990,16 +2007,20 @@ algorithm
     end if;
 
     // translate adjacency matrix to BackendDAE sparsity pattern
-    intsp := BackendDAEUtil.absIncidenceMatrix(Util.getOption(system.m));
+    intsp := BackendDAEUtil.incidenceMatrixToSparsePattern(Util.getOption(system.m));
     intspList := arrayList(intsp);
-    intspT := BackendDAEUtil.absIncidenceMatrix(Util.getOption(system.mT));
+    intspT := BackendDAEUtil.incidenceMatrixToSparsePattern(Util.getOption(system.mT));
     intspTList := arrayList(intspT);
-
     // get number of non zero elements
     nonZeroElements := List.lengthListElements(intspList);
+    if debug then
+      BackendDump.dumpSparsePatternArray(intsp);
+      print("daeMode prepared a pattern: " + realString(clock()) + "\n");
+      BackendDump.dumpSparsePatternArray(intspT);
+      print("daeMode prepared a pattern: " + realString(clock()) + "\n");
+    end if;
 
     // translated to DAE.ComRefs
-    varCrefsList := List.map(BackendVariable.varList(inVars), BackendVariable.varCref);
     varCrefs := listArray(varCrefsList);
 
     translated := list(list(arrayGet(varCrefs, i) for i in lst) for lst in intspList);
@@ -2010,14 +2031,14 @@ algorithm
     // build sparse pattern
     bdaeSP :=(sparsetuple, sparsetupleT, (varCrefsList, varCrefsList), nonZeroElements);
 
+    // get coloring
     coloredArray := SymbolicJacobian.createColoring(intsp, intspT, listLength(varCrefsList), listLength(varCrefsList));
     coloring := list(list(arrayGet(varCrefs, i) for i in lst) for lst in coloredArray);
-    // get coloring
     // or without coloring
-    //coloring = List.transposeList({inDepCompRefs});
+    //coloring := List.transposeList({varCrefsList});
 
     // translate to SimCode sparsity
-    ({simCodeJac}, _) := createSymbolicJacobianssSimCode({(NONE(), bdaeSP, coloring)}, crefToSimVarHT, 0, {"daeMode"});
+    (simCodeJac, _) :=  createSimCodeSparsePattenDAEmode({(NONE(), bdaeSP, coloring)}, crefToSimVarHT, 0, "daeMode");
 
     outSimCodeJac := SOME(simCodeJac);
   else
@@ -2039,6 +2060,20 @@ algorithm
   end if;
 end getSimCodeDAEModeDataEqns;
 
+protected function sortSimVarsAndWriteIndex
+  input list<SimCodeVar.SimVar> inSimVar;
+  input SimCode.HashTableCrefToSimVar crefToSimVarHT;
+  output list<SimCodeVar.SimVar> outSimVar;
+protected
+  list<DAE.ComponentRef> crefs;
+  list<SimCodeVar.SimVar> sorted;
+algorithm
+  crefs := List.map(inSimVar, getSimVarCompRef);
+  sorted := getSimVars2Crefs(crefs, crefToSimVarHT);
+  sorted := List.sort(sorted, compareVarIndexGt);
+  outSimVar := rewriteIndex(sorted, 0);
+end sortSimVarsAndWriteIndex;
+
 // =============================================================================
 // section for zeroCrossingsEquations
 //
@@ -4526,6 +4561,94 @@ algorithm
   end if;
 end sortSparsePattern;
 
+protected function sortSparsePatternDAEmode
+  input list<SimCodeVar.SimVar> inSimVars;
+  input list<tuple<DAE.ComponentRef, list<DAE.ComponentRef>>> inSparsePattern;
+  input Boolean useFMIIndex;
+  output list<tuple<Integer, list<Integer>>> outSparse = {};
+protected
+  HashTable.HashTable ht;
+  DAE.ComponentRef cref;
+  Integer size, i, j;
+  list<Integer> intLst;
+  list<DAE.ComponentRef> crefs;
+algorithm
+  //create HT
+  size := listLength(inSimVars);
+  if size>0 then
+    ht := HashTable.emptyHashTableSized(size);
+    for var in inSimVars loop
+      if not useFMIIndex then
+        SimCodeVar.SIMVAR(name = cref, index=i) := var;
+      else
+        SimCodeVar.SIMVAR(name = cref) := var;
+        i := getVariableIndex(var);
+      end if;
+      //print("Setup HashTable with cref: " + ComponentReference.printComponentRefStr(cref) + " index: "+ intString(i) + "\n");
+      ht := BaseHashTable.add((cref, i), ht);
+    end for;
+
+    //translate
+    i := 0;
+    for tpl in inSparsePattern loop
+       (cref, crefs) := tpl;
+       intLst := {};
+       for cr in crefs loop
+         j := BaseHashTable.get(cr, ht);
+         intLst := j :: intLst;
+       end for;
+       intLst := List.sort(intLst, intGt);
+       outSparse := (i, intLst) :: outSparse;
+       i := i + 1;
+    end for;
+    outSparse := listReverse(outSparse);
+  end if;
+end sortSparsePatternDAEmode;
+
+protected function sortSparsePatternTDAEmode
+  input list<SimCodeVar.SimVar> inSimVars;
+  input list<tuple<DAE.ComponentRef, list<DAE.ComponentRef>>> inSparsePattern;
+  input Boolean useFMIIndex;
+  output list<tuple<Integer, list<Integer>>> outSparse = {};
+protected
+  HashTable.HashTable ht;
+  DAE.ComponentRef cref;
+  Integer size, i, j;
+  list<Integer> intLst;
+  list<DAE.ComponentRef> crefs;
+algorithm
+  //create HT
+  size := listLength(inSimVars);
+  if size>0 then
+    ht := HashTable.emptyHashTableSized(size);
+    for var in inSimVars loop
+      if not useFMIIndex then
+        SimCodeVar.SIMVAR(name = cref, index=i) := var;
+      else
+        SimCodeVar.SIMVAR(name = cref) := var;
+        i := getVariableIndex(var);
+      end if;
+      //print("Setup HashTable with cref: " + ComponentReference.printComponentRefStr(cref) + " index: "+ intString(i) + "\n");
+      ht := BaseHashTable.add((cref, i), ht);
+    end for;
+
+    //translate
+    for tpl in inSparsePattern loop
+       (cref, crefs) := tpl;
+       i := BaseHashTable.get(cref, ht);
+       intLst := {};
+       j := 0;
+       for cr in crefs loop
+         intLst := j :: intLst;
+         j := j + 1;
+       end for;
+       intLst := listReverse(intLst);
+       outSparse := (i, intLst) :: outSparse;
+    end for;
+  outSparse := List.sort(outSparse, Util.compareTupleIntGt);
+  end if;
+end sortSparsePatternTDAEmode;
+
 protected function sortColoring
   input list<SimCodeVar.SimVar> inSimVars;
   input list<list<DAE.ComponentRef>> inColoring;
@@ -4584,6 +4707,103 @@ algorithm
   end for;
 end dumpSparsePattern;
 
+protected function createSimCodeSparsePattenDAEmode
+"fuction translates the sparse pattern of the daeMode
+ author: wbraun"
+  input BackendDAE.SymbolicJacobians inSymJacobian;
+  input SimCode.HashTableCrefToSimVar inSimVarHT;
+  input Integer iuniqueEqIndex;
+  input String inName;
+  output SimCode.JacobianMatrix outJacobianMatrix;
+  output Integer ouniqueEqIndex;
+algorithm
+  (outJacobianMatrix, ouniqueEqIndex) :=
+  matchcontinue (inSymJacobian)
+    local
+      list<DAE.ComponentRef> diffCompRefs, diffedCompRefs;
+
+      Integer uniqueEqIndex;
+
+      String name, s, dummyVar;
+
+      list<SimCodeVar.SimVar> seedVars, indexVars, seedIndexVars;
+
+      list<tuple<DAE.ComponentRef, list<DAE.ComponentRef>>> sparsepattern, sparsepatternT;
+      list<list<DAE.ComponentRef>> colsColors;
+      Integer maxColor;
+
+      list<tuple<Integer, list<Integer>>> sparseInts, sparseIntsT;
+      array<tuple<Integer, list<Integer>>> sparseArrayT;
+      list<list<Integer>> coloring;
+      Option<BackendDAE.SymbolicJacobian> optionBDAE;
+
+      constant Boolean debug = false;
+
+    // if only sparsity pattern is generated
+    case ({(optionBDAE, (sparsepattern, sparsepatternT, (diffCompRefs, diffedCompRefs), _), colsColors)})
+      guard  checkForEmptyBDAE(optionBDAE)
+      equation
+        if debug then
+          print("Start sparse pattern without analytical Jacobians\n");
+        end if;
+
+        seedVars = getSimVars2Crefs(diffCompRefs, inSimVarHT);
+        seedVars = List.sort(seedVars, compareVarIndexGt);
+
+        if debug then
+          print("diffCrefs: " + ComponentReference.printComponentRefListStr(diffCompRefs) + "\n");
+          print("\n---+++  seedVars +++---\n");
+          print(Tpl.tplString(SimCodeDump.dumpVarsShort, seedVars));
+        end if;
+
+        indexVars = getSimVars2Crefs(diffedCompRefs, inSimVarHT);
+        indexVars = List.sort(indexVars, compareVarIndexGt);
+
+        if debug then
+          print("diffedCrefs: " + ComponentReference.printComponentRefListStr(diffedCompRefs) + "\n");
+          print("\n---+++  indexVars +++---\n");
+          print(Tpl.tplString(SimCodeDump.dumpVarsShort, indexVars));
+          print("\n---+++  sparse pattern vars +++---\n");
+          dumpSparsePattern(sparsepattern);
+          print("\n---+++  sparse pattern transpose +++---\n");
+          dumpSparsePattern(sparsepatternT);
+        end if;
+        seedVars = rewriteIndex(seedVars, 0);
+        indexVars = rewriteIndex(indexVars, 0);
+        seedIndexVars = listAppend(seedVars, indexVars);
+
+        //sparseInts = sortSparsePattern(seedIndexVars, sparsepattern, false);
+        //sparseIntsT = sortSparsePattern(seedIndexVars, sparsepatternT, false);
+        sparseInts = sortSparsePatternDAEmode(seedIndexVars, sparsepattern, false);
+
+        sparseArrayT = arrayCreate(listLength(sparseInts), (-1,{}));
+        sparseArrayT = SymbolicJacobian.transposeSparsePatternTuple(sparseInts, sparseArrayT);
+        sparseIntsT = arrayList(sparseArrayT);
+
+        maxColor = listLength(colsColors);
+        s = intString(listLength(diffedCompRefs));
+        coloring = sortColoring(seedVars, colsColors);
+
+        if debug then
+          print("analytical Jacobians -> transformed to SimCode for Matrix " + inName + " time: " + realString(clock()) + "\n");
+          print("\n---+++  sparse pattern vars +++---\n");
+          dumpSparsePatternInt(sparseInts);
+          print("\n---+++  sparse pattern transpose +++---\n");
+          dumpSparsePatternInt(sparseIntsT);
+        end if;
+
+     then
+        (({(({}, {}, s))}, {}, inName, (sparseInts, sparseIntsT), coloring, maxColor, -1), iuniqueEqIndex);
+
+    else
+      equation
+        Error.addInternalError("Function createSimCodeSparsePattenDAEmode failed", sourceInfo());
+      then
+        fail();
+  end matchcontinue;
+end createSimCodeSparsePattenDAEmode;
+
+
 // =============================================================================
 // section with unsorted function
 //
@@ -9146,6 +9366,18 @@ algorithm
   result := index1 > index2;
 end compareVarIndexGt;
 
+public function compareSimVarTupleIndexGt
+  input tuple<SimCodeVar.SimVar, BackendDAE.Var> var1;
+  input tuple<SimCodeVar.SimVar, BackendDAE.Var> var2;
+  output Boolean result;
+protected
+  Integer index1, index2;
+algorithm
+  (SimCodeVar.SIMVAR(variable_index=SOME(index1)),_) := var1;
+  (SimCodeVar.SIMVAR(variable_index=SOME(index2)),_) := var2;
+  result := index1 > index2;
+end compareSimVarTupleIndexGt;
+
 public function countDynamicExternalFunctions
   input list<SimCode.Function> inFncLst;
   output Integer outDynLoadFuncs;
diff --git a/Compiler/Template/CodegenC.tpl b/Compiler/Template/CodegenC.tpl
index a2d1d006bb4..b2a2add9a2b 100644
--- a/Compiler/Template/CodegenC.tpl
+++ b/Compiler/Template/CodegenC.tpl
@@ -836,6 +836,13 @@ template simulationFile_dae_header(SimCode simCode)
     ;separator="\n"%>
     >>
     /* adrpo: leave a newline at the end of file to get rid of the warning */
+    case simCode as SIMCODE(__) then
+    <<
+    #ifndef <%fileNamePrefix%>_16DAE_H
+    #define <%fileNamePrefix%>_16DAE_H
+    #endif
+    <%\n%>
+    >>
   end match
 end simulationFile_dae_header;
 
@@ -847,9 +854,9 @@ template simulationFile_dae(SimCode simCode)
         daeModeData=SOME(DAEMODEDATA(daeEquations=daeEquations, sparsityPattern=sparsityPattern,
                                      algebraicDAEVars=algebraicDAEVars, residualVars=residualVars))) then
      let modelNamePrefixStr = modelNamePrefix(simCode)
-     let initDAEmode = 
+     let initDAEmode =
        match sparsityPattern
-       case SOME((_, _, _, (sparse,_), colorList, maxColor, _)) then
+       case SOME((_, _, _, (_, sparse), colorList, maxColor, _)) then
          '<%initializeDAEmodeData(listLength(residualVars), listLength(algebraicDAEVars), sparse, colorList, maxColor, modelNamePrefixStr)%>'
        case NONE() then
          'int <%symbolName(modelNamePrefixStr,"initializeDAEmodeData")%>(DATA *inData, DAEMODE_DATA* daeModeData){ return -1; }'
@@ -3807,7 +3814,7 @@ template evaluateDAEResiduals(list<list<SimEqSystem>> resEquations, String model
 
   <%systems%>
   /* for residuals DAE variables */
-  int <%symbolName(modelNamePrefix,"evaluateDAEResiduals")%>(DATA *data, threadData_t *threadData, double* residualVars)
+  int <%symbolName(modelNamePrefix,"evaluateDAEResiduals")%>(DATA *data, threadData_t *threadData)
   {
     TRACE_PUSH
     data->simulationInfo->callStatistics.functionEvalDAE++;
diff --git a/SimulationRuntime/c/simulation/solver/ida_solver.c b/SimulationRuntime/c/simulation/solver/ida_solver.c
index 0ec81e4757d..7991639a39d 100644
--- a/SimulationRuntime/c/simulation/solver/ida_solver.c
+++ b/SimulationRuntime/c/simulation/solver/ida_solver.c
@@ -81,9 +81,7 @@ static int jacobiancoloredKLUNum(realtype tt, realtype cj,
     N_Vector yy, N_Vector yp, N_Vector rr, SlsMat Jac, void *user_data,
     N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
 
-
 static int residualFunctionIDA(double time, N_Vector yy, N_Vector yp, N_Vector res, void* userData);
-static int residualFunctionIDADAEmode(double time, N_Vector yy, N_Vector yp, N_Vector res, void* userData);
 int rootsFunctionIDA(double time, N_Vector yy, N_Vector yp, double *gout, void* userData);
 
 int checkIDAflag(int flag)
@@ -149,7 +147,6 @@ ida_solver_initial(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo
     if (compiledInDAEMode)
     {
       idaData->daeMode = 1;
-      idaData->residualFunction = residualFunctionIDADAEmode;
       idaData->N = data->modelData->nStates + data->simulationInfo->daeModeData->nAlgebraicDAEVars;
     }
     else
@@ -260,7 +257,7 @@ ida_solver_initial(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo
   }
   else
   {
-    idaData->linearSolverMethod = IDA_LS_DENSE;
+    idaData->linearSolverMethod = IDA_LS_KLU;
   }
 
   switch (idaData->linearSolverMethod){
@@ -317,31 +314,22 @@ ida_solver_initial(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo
   }
   else
   {
-    /* in dae mode use for now internal Jacobian calculation */
-    if (idaData->daeMode)
+    if (idaData->linearSolverMethod == IDA_LS_KLU)
     {
-      idaData->jacobianMethod = NUMJAC;
-      warningStreamPrint(LOG_SOLVER, 0, "Running in DAE mode so currently no colored jacobian available, yet!");
+      idaData->jacobianMethod = KLUSPARSE;
     }
     else
     {
-      if (idaData->linearSolverMethod == IDA_LS_KLU)
-      {
-        idaData->jacobianMethod = KLUSPARSE;
-      }
-      else
-      {
-        idaData->jacobianMethod = COLOREDNUMJAC;
-      }
+      idaData->jacobianMethod = COLOREDNUMJAC;
     }
   }
 
   /* selects the calculation method of the jacobian */
-  if(idaData->jacobianMethod == COLOREDNUMJAC ||
+/* in daeMode sparse pattern is already initialized in DAEMODE_DATA */
+  if(!idaData->daeMode && (idaData->jacobianMethod == COLOREDNUMJAC ||
       idaData->jacobianMethod == COLOREDSYMJAC ||
       idaData->jacobianMethod == KLUSPARSE ||
-      idaData->jacobianMethod == KLUCOLORED ||
-      idaData->jacobianMethod == SYMJAC)
+      idaData->jacobianMethod == SYMJAC))
   {
     if (data->callback->initialAnalyticJacobianA(data, threadData))
     {
@@ -353,18 +341,28 @@ ida_solver_initial(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo
   /* set up the appropriate function pointer */
   if (idaData->linearSolverMethod == IDA_LS_KLU)
   {
-    flag = IDAKLU(idaData->ida_mem, idaData->N, data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A].sparsePattern.numberOfNoneZeros);
-    idaData->tmpJac = NewSparseMat(data->modelData->nStates, data->modelData->nStates, data->modelData->nStates);
+    if (idaData->daeMode)
+    {
+      idaData->NNZ = data->simulationInfo->daeModeData->sparsePattern->numberOfNoneZeros;
+      flag = IDAKLU(idaData->ida_mem, idaData->N, idaData->NNZ);
+    }
+    else
+    {
+      idaData->NNZ = data->simulationInfo->analyticJacobians[data->callback->INDEX_JAC_A].sparsePattern.numberOfNoneZeros;
+      flag = IDAKLU(idaData->ida_mem, idaData->N, idaData->NNZ);
+      /* to add a cj identety matrix */
+      idaData->tmpJac = NewSparseMat(idaData->N, idaData->N, idaData->NNZ);
+    }
+    if (checkIDAflag(flag)){
+      throwStreamPrint(threadData, "##IDA## Setup of linear solver KLU failed!");
+    }
+
     switch (idaData->jacobianMethod){
     case KLUSPARSE:
       flag = IDASlsSetSparseJacFn(idaData->ida_mem, jacobianSparseNum);
       break;
-    case KLUCOLORED:
-      idaData->denseJac = NewDenseMat(data->modelData->nStates, data->modelData->nStates);
-      flag = IDASlsSetSparseJacFn(idaData->ida_mem, jacobiancoloredKLUNum);
-      break;
     default:
-      throwStreamPrint(threadData,"unrecognized jacobian calculation method %s", (const char*)omc_flagValue[FLAG_JACOBIAN]);
+      throwStreamPrint(threadData,"For the klu solver jacobian calculation method has to be %s", JACOBIAN_METHOD[KLUSPARSE]);
       break;
     }
   }
@@ -482,13 +480,9 @@ ida_solver_deinitial(IDA_SOLVER *idaData){
   free(idaData->ysave);
   free(idaData->ypsave);
   free(idaData->delta_hh);
-  if (idaData->linearSolverMethod == IDA_LS_KLU){
+  if (!idaData->daeMode && idaData->linearSolverMethod == IDA_LS_KLU){
     DestroySparseMat(idaData->tmpJac);
   }
-  if (idaData->jacobianMethod == KLUCOLORED)
-  {
-    DestroyMat(idaData->denseJac);
-  }
 
   if (idaData->daeMode)
   {
@@ -525,7 +519,7 @@ ida_solver_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo)
 
   IDA_SOLVER *idaData = (IDA_SOLVER*) solverInfo->solverData;
 
-  SIMULATION_DATA *sData = data->localData[0];
+    SIMULATION_DATA *sData = data->localData[0];
   SIMULATION_DATA *sDataOld = data->localData[1];
   MODEL_DATA *mData = (MODEL_DATA*) data->modelData;
 
@@ -545,6 +539,8 @@ ida_solver_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo)
   /* reinit solver */
   if (!idaData->setInitialSolution)
   {
+    debugStreamPrint(LOG_SOLVER, 0, "Re-initialized IDA Solver");
+
     /* initialize states and der(states) */
     if (idaData->daeMode)
     {
@@ -552,14 +548,13 @@ ida_solver_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo)
       /* and  also algebraic vars */
       data->simulationInfo->daeModeData->getAlgebraicDAEVars(data, threadData, idaData->states + data->modelData->nStates);
       memcpy(idaData->statesDer, data->localData[1]->realVars + data->modelData->nStates, sizeof(double)*data->modelData->nStates);
+
     }
     flag = IDAReInit(idaData->ida_mem,
         solverInfo->currentTime,
         idaData->y,
         idaData->yp);
 
-    debugStreamPrint(LOG_SOLVER, 0, "Re-initialized IDA Solver");
-
     if (checkIDAflag(flag)){
       throwStreamPrint(threadData, "##IDA## Something goes wrong while reinit IDA solver after event!");
     }
@@ -769,122 +764,68 @@ int residualFunctionIDA(double time, N_Vector yy, N_Vector yp, N_Vector res, voi
   {
     setContext(data, &time, CONTEXT_ODE);
   }
-  printCurrentStatesVector(LOG_DASSL_STATES, states, data, time);
-  printVector(LOG_DASSL_STATES, "yd", statesDer, data->modelData->nStates, time);
-
   timeBackup = data->localData[0]->timeValue;
   data->localData[0]->timeValue = time;
 
   saveJumpState = threadData->currentErrorStage;
   threadData->currentErrorStage = ERROR_INTEGRATOR;
 
-
   /* try */
 #if !defined(OMC_EMCC)
   MMC_TRY_INTERNAL(simulationJumpBuffer)
 #endif
-  data->simulationInfo->daeModeData->setAlgebraicDAEVars(data, threadData, idaData->states + data->modelData->nStates);
+
   /* if sensitivity mode update also bound parameters*/
   if (idaData->idaSmode)
   {
     data->callback->updateBoundParameters(data, threadData);
   }
-
-  /* read input vars */
-  externalInputUpdate(data);
-  data->callback->input_function(data, threadData);
-
-  /* eval input vars */
-  data->callback->functionODE(data, threadData);
-
-  /* get the difference between the temp_xd(=localData->statesDerivatives)
-     and xd(=statesDerivativesBackup) */
-  for(i=0; i < data->modelData->nStates; i++)
-  {
-    NV_Ith_S(res, i) = data->localData[0]->realVars[data->modelData->nStates + i] - NV_Ith_S(yp, i);
-  }
-  printVector(LOG_DASSL_STATES, "dd", delta, data->modelData->nStates, time);
-  success = 1;
-#if !defined(OMC_EMCC)
-  MMC_CATCH_INTERNAL(simulationJumpBuffer)
-#endif
-
-  if (!success) {
-    retVal = -1;
-  }
-
-  threadData->currentErrorStage = saveJumpState;
-
-  data->localData[0]->timeValue = timeBackup;
-
-  if (data->simulationInfo->currentContext == CONTEXT_ODE){
-    unsetContext(data);
-  }
-
-  TRACE_POP
-  return retVal;
-}
-
-
-int residualFunctionIDADAEmode(double time, N_Vector yy, N_Vector yp, N_Vector res, void* userData)
-{
-  TRACE_PUSH
-  IDA_SOLVER* idaData = (IDA_SOLVER*)userData;
-  DATA* data = (DATA*)(((IDA_USERDATA*)idaData->simData)->data);
-  threadData_t* threadData = (threadData_t*)(((IDA_USERDATA*)((IDA_SOLVER*)userData)->simData)->threadData);
-
-  double timeBackup;
-  long int i;
-  int saveJumpState;
-  int success = 0, retVal = 0;
-  double *states = N_VGetArrayPointer(yy);
-  double *statesDer = N_VGetArrayPointer(yp);
-  double *delta  = N_VGetArrayPointer(res);
-
-  /* context */
-  if (data->simulationInfo->currentContext == CONTEXT_ALGEBRAIC)
+  /* if daeMode update also all dynamic algebraic equations */
+  if (idaData->daeMode)
   {
-    setContext(data, &time, CONTEXT_ODE);
+    /* set state, state derivative and dynamic algebraic
+    * variables for evaluateDAEResiduals evaluation
+    */
+    memcpy(data->localData[0]->realVars, states, sizeof(double)*data->modelData->nStates);
+    memcpy(data->localData[0]->realVars + data->modelData->nStates, statesDer, sizeof(double)*data->modelData->nStates);
+    data->simulationInfo->daeModeData->setAlgebraicDAEVars(data, threadData, states + data->modelData->nStates);
   }
 
   /* debug */
   if (ACTIVE_STREAM(LOG_DASSL_STATES)){
-    printCurrentStatesVector(LOG_DASSL_STATES, states, data, time);
-    printVector(LOG_DASSL_STATES, "yprime", statesDer, data->modelData->nStates, time);
+    printCurrentStatesVector(LOG_DASSL_STATES, data->localData[0]->realVars, data, time);
+    printVector(LOG_DASSL_STATES, "yprime", data->localData[0]->realVars + data->modelData->nStates, data->modelData->nStates, time);
+    if (idaData->daeMode)
+    {
+      printVector(LOG_DASSL_STATES, "yalg", states + data->modelData->nStates, data->simulationInfo->daeModeData->nAlgebraicDAEVars, time);
+    }
   }
 
-  /* set time */
-  timeBackup = data->localData[0]->timeValue;
-  data->localData[0]->timeValue = time;
-
-  /* set state, state derivative and dynamic algebraic
-   * variables for evaluateDAEResiduals evaluation
-   */
-  memcpy(data->localData[0]->realVars, idaData->states, sizeof(double)*data->modelData->nStates);
-  memcpy(data->localData[0]->realVars + data->modelData->nStates, idaData->statesDer, sizeof(double)*data->modelData->nStates);
-  data->simulationInfo->daeModeData->setAlgebraicDAEVars(data, threadData, idaData->states + data->modelData->nStates);
-
-  saveJumpState = threadData->currentErrorStage;
-  threadData->currentErrorStage = ERROR_INTEGRATOR;
-
-  /* try */
-#if !defined(OMC_EMCC)
-  MMC_TRY_INTERNAL(simulationJumpBuffer)
-#endif
-
   /* read input vars */
   externalInputUpdate(data);
   data->callback->input_function(data, threadData);
 
-  /* eval residual vars */
-  data->simulationInfo->daeModeData->evaluateDAEResiduals(data, threadData);
-
-  /* get data->simulationInfo->residualVars  */
-  for(i=0; i < idaData->N; i++)
+  if (idaData->daeMode)
   {
-    NV_Ith_S(res, i) = data->simulationInfo->daeModeData->residualVars[i];
+    /* eval residual vars */
+    data->simulationInfo->daeModeData->evaluateDAEResiduals(data, threadData);
+    /* get residual variables */
+    for(i=0; i < idaData->N; i++)
+    {
+      NV_Ith_S(res, i) = data->simulationInfo->daeModeData->residualVars[i];
+    }
   }
-  printVector(LOG_DASSL_STATES, "residual", delta, idaData->N, time);
+  else
+  {
+    /* eval function ODE */
+    data->callback->functionODE(data, threadData);
+    for(i=0; i < idaData->N; i++)
+    {
+        NV_Ith_S(res, i) = data->localData[0]->realVars[data->modelData->nStates + i] - NV_Ith_S(yp, i);
+    }
+  }
+
+  printVector(LOG_DASSL_STATES, "delta", delta, idaData->N, time);
   success = 1;
 #if !defined(OMC_EMCC)
   MMC_CATCH_INTERNAL(simulationJumpBuffer)
@@ -957,7 +898,8 @@ int rootsFunctionIDA(double time, N_Vector yy, N_Vector yp, double *gout, void*
 
 /*
  *  function calculates a jacobian matrix by
- *  numerical method finite differences
+ *  numerical method finite differences with coloring
+ *  into a dense DlsMat matrix
  */
 static
 int jacOwnNumColoredIDA(double tt, N_Vector yy, N_Vector yp, N_Vector rr, DlsMat Jac, double cj, void *userData)
@@ -989,19 +931,30 @@ int jacOwnNumColoredIDA(double tt, N_Vector yy, N_Vector yp, N_Vector rr, DlsMat
   IDAGetCurrentStep(ida_mem, &currentStep);
   IDAGetErrWeights(ida_mem, idaData->errwgt);
 
+  SPARSE_PATTERN* sparsePattern;
+
+  /* set sparse pattern */
+  if (idaData->daeMode)
+  {
+    sparsePattern = data->simulationInfo->daeModeData->sparsePattern;
+  }
+  else
+  {
+    sparsePattern = &(data->simulationInfo->analyticJacobians[index].sparsePattern);
+  }
+
   setContext(data, &tt, CONTEXT_JACOBIAN);
 
-  for(i = 0; i < data->simulationInfo->analyticJacobians[index].sparsePattern.maxColors; i++)
+  for(i = 0; i < sparsePattern->maxColors; i++)
   {
-    for(ii=0; ii < data->simulationInfo->analyticJacobians[index].sizeCols; ii++)
+    for(ii=0; ii < idaData->N; ii++)
     {
-      if(data->simulationInfo->analyticJacobians[index].sparsePattern.colorCols[ii]-1 == i)
+      if(sparsePattern->colorCols[ii]-1 == i)
       {
         delta_hhh = currentStep * yprime[ii];
         delta_hh[ii] = delta_h * fmax(fmax(fabs(states[ii]),fabs(delta_hhh)),fabs(1./errwgt[ii]));
         delta_hh[ii] = (delta_hhh >= 0 ? delta_hh[ii] : -delta_hh[ii]);
         delta_hh[ii] = (states[ii] + delta_hh[ii]) - states[ii];
-
         ysave[ii] = states[ii];
         states[ii] += delta_hh[ii];
 
@@ -1018,20 +971,33 @@ int jacOwnNumColoredIDA(double tt, N_Vector yy, N_Vector yp, N_Vector rr, DlsMat
 
     increaseJacContext(data);
 
-    for(ii = 0; ii < data->simulationInfo->analyticJacobians[index].sizeCols; ii++)
+    for(ii = 0; ii < idaData->N; ii++)
     {
-      if(data->simulationInfo->analyticJacobians[index].sparsePattern.colorCols[ii]-1 == i)
+      if(sparsePattern->colorCols[ii]-1 == i)
       {
-        if(ii==0)
-          j = 0;
+        if (idaData->daeMode)
+        {
+          j = sparsePattern->leadindex[ii];
+          while(j < sparsePattern->leadindex[ii+1])
+          {
+            l  =  sparsePattern->index[j];
+            DENSE_ELEM(Jac, l, ii) = (newdelta[l] - delta[l]) * delta_hh[ii];
+            j++;
+          };
+        }
         else
-          j = data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[ii-1];
-        while(j < data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[ii])
         {
-          l  =  data->simulationInfo->analyticJacobians[index].sparsePattern.index[j];
-          DENSE_ELEM(Jac, l, ii) = (newdelta[l] - delta[l]) * delta_hh[ii];
-          j++;
-        };
+          if(ii==0)
+            j = 0;
+          else
+            j = sparsePattern->leadindex[ii-1];
+          while(j < sparsePattern->leadindex[ii])
+          {
+            l  =  sparsePattern->index[j];
+            DENSE_ELEM(Jac, l, ii) = (newdelta[l] - delta[l]) * delta_hh[ii];
+            j++;
+          };
+        }
         states[ii] = ysave[ii];
         if (idaData->daeMode)
         {
@@ -1047,7 +1013,10 @@ int jacOwnNumColoredIDA(double tt, N_Vector yy, N_Vector yp, N_Vector rr, DlsMat
 }
 
 /*
- * provides a numerical Jacobian to be used with DASSL
+ * Wrapper function jacOwnNumColoredIDA
+ *  function calculates a jacobian matrix by
+ *  numerical method finite differences with coloring
+ *  into a dense DlsMat matrix
  */
 static int jacobianOwnNumColoredIDA(long int Neq, double tt, double cj,
     N_Vector yy, N_Vector yp, N_Vector rr,
@@ -1088,7 +1057,8 @@ static int jacobianOwnNumColoredIDA(long int Neq, double tt, double cj,
 
 /*
  *  function calculates a jacobian matrix by
- *  numerical method finite differences
+ *  numerical method finite differences without coloring
+ *  into a dense DlsMat matrix
  */
 static
 int jacOwnNumIDA(double tt, N_Vector yy, N_Vector yp, N_Vector rr, DlsMat Jac, double cj, void *userData)
@@ -1142,6 +1112,7 @@ int jacOwnNumIDA(double tt, N_Vector yy, N_Vector yp, N_Vector rr, DlsMat Jac, d
     for(j = 0; j < idaData->N; j++)
     {
       DENSE_ELEM(Jac, j, i) = (newdelta[j] - delta[j]) * deltaInv;
+
     }
     states[i] = ysave;
     if (idaData->daeMode){
@@ -1194,54 +1165,7 @@ static int jacobianOwnNumIDA(long int Neq, double tt, double cj,
   return 0;
 }
 
-static void transposeJac(SlsMat spJac)
-{
-  int i,j,k,index;
-  const int N = spJac->N;
-  const int NNZ = spJac->NNZ;
-  int *count = (int*) calloc(N, sizeof(int));
-
-  /*  Initialized to zero. */
-  SlsMat tmpJac = NewSparseMat(N, N, NNZ);
-
-  /* First find the column lengths for spJac^{T}
-   * i.e. the row lengths of A.
-   * Temporary counters for each row of A.
-   */
-  for (i=0; i<N; i++)
-  {
-    for (j=spJac->colptrs[i];j<spJac->colptrs[i+1];j++)
-    {
-      k=spJac->rowvals[j];
-      count[k]++;
-    }
-  }
-
-  /* Now set spJac->colptrs. 0th entry stays 0. */
-  tmpJac->colptrs[0] = 0;
-  for (j=0;j<N;j++)
-  {
-    tmpJac->colptrs[j+1]=tmpJac->colptrs[j]+count[j];
-    count[j]=0;
-  }
-  /* Main loop.*/
-  for (i=0;i<N;i++)
-  {
-    for (j=spJac->colptrs[i];j<spJac->colptrs[i+1];j++)
-    {
-      k=spJac->rowvals[j];
-      /*Element’s position in column of Jac^T .*/
-      index=tmpJac->colptrs[k]+count[k];
-      tmpJac->rowvals[index]=i;
-      tmpJac->data[index]=spJac->data[j];
-      /*Increment counter for next element in that column. */
-      count[k]++;
-    }
-  }
-  CopySparseMat(tmpJac,spJac);
-  free(count);
-}
-
+/* Element function for sparse matrix set */
 static void setJacElementKluSparse(int row, int col, double value, int nth, SlsMat spJac)
 {
   if (col > 0){
@@ -1255,7 +1179,8 @@ static void setJacElementKluSparse(int row, int col, double value, int nth, SlsM
 
 /*
  *  function calculates a jacobian matrix by
- *  numerical method finite differences
+ *  numerical method finite differences with coloring
+ *  into a sparse SlsMat matrix
  */
 static
 int jacobianSparseNumIDA(double tt, N_Vector yy, N_Vector yp, N_Vector rr, SlsMat Jac, double cj, void *userData)
@@ -1273,10 +1198,13 @@ int jacobianSparseNumIDA(double tt, N_Vector yy, N_Vector yp, N_Vector rr, SlsMa
   double *newdelta = N_VGetArrayPointer(idaData->newdelta);
   double *errwgt = N_VGetArrayPointer(idaData->errwgt);
 
-  double ysave, ypsave;
+  SPARSE_PATTERN* sparsePattern;
+
+  double *ysave = idaData->ysave;
+  double *ypsave = idaData->ypsave;
 
   double delta_h = idaData->sqrteps;
-  double delta_hh;
+  double *delta_hh = idaData->delta_hh;
   double delta_hhh;
   double deltaInv;
 
@@ -1289,50 +1217,81 @@ int jacobianSparseNumIDA(double tt, N_Vector yy, N_Vector yp, N_Vector rr, SlsMa
   IDAGetCurrentStep(ida_mem, &currentStep);
   IDAGetErrWeights(ida_mem, idaData->errwgt);
 
+  /* set sparse pattern */
+  if (idaData->daeMode)
+  {
+    sparsePattern = data->simulationInfo->daeModeData->sparsePattern;
+  }
+  else
+  {
+    sparsePattern = &(data->simulationInfo->analyticJacobians[index].sparsePattern);
+  }
+
   /* it's needed to clear the matrix */
   SlsSetToZero(Jac);
 
   setContext(data, &tt, CONTEXT_JACOBIAN);
 
-  for(i = 0; i < idaData->N; i++)
+  for(i = 0; i < sparsePattern->maxColors; i++)
   {
-    delta_hhh = currentStep * yprime[i];
-    delta_hh = delta_h * fmax(fmax(fabs(states[i]),fabs(delta_hhh)),fabs(1./errwgt[i]));
-    delta_hh = (delta_hhh >= 0 ? delta_hh : -delta_hh);
-    delta_hh = (states[i] + delta_hh) - states[i];
-    deltaInv = 1. / delta_hh;
-    ysave = states[i];
-    states[i] += delta_hh;
-    if (idaData->daeMode){
-      ypsave = yprime[i];
-      yprime[i] += cj * delta_hh;
+    for(ii=0; ii < idaData->N; ii++)
+    {
+      if(sparsePattern->colorCols[ii]-1 == i)
+      {
+        delta_hhh = currentStep * yprime[ii];
+        delta_hh[ii] = delta_h * fmax(fmax(fabs(states[ii]),fabs(delta_hhh)),fabs(1./errwgt[ii]));
+        delta_hh[ii] = (delta_hhh >= 0 ? delta_hh[ii] : -delta_hh[ii]);
+        delta_hh[ii] = (states[ii] + delta_hh[ii]) - states[ii];
+        ysave[ii] = states[ii];
+        states[ii] += delta_hh[ii];
+
+        if (idaData->daeMode){
+          ypsave[ii] = yprime[ii];
+          yprime[ii] += cj * delta_hh[ii];
+        }
+
+        delta_hh[ii] = 1. / delta_hh[ii];
+      }
     }
 
     (*idaData->residualFunction)(tt, yy, yp, idaData->newdelta, userData);
 
     increaseJacContext(data);
 
-    ii = (i == 0) ?  0 : data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[i-1];
-
-    while(ii < data->simulationInfo->analyticJacobians[index].sparsePattern.leadindex[i])
+    for(ii = 0; ii < idaData->N; ii++)
     {
-      j  =  data->simulationInfo->analyticJacobians[index].sparsePattern.index[ii];
-      setJacElementKluSparse(j, i, (newdelta[j] - delta[j]) * deltaInv, nth, Jac);
-      ii++;
-      nth++;
-    };
-
-    states[i] = ysave;
-    if (idaData->daeMode)
-    {
-      yprime[i] = ypsave;
+      if(sparsePattern->colorCols[ii]-1 == i)
+      {
+        if (idaData->daeMode)
+        {
+          nth = sparsePattern->leadindex[ii];
+          while(nth < sparsePattern->leadindex[ii+1])
+          {
+            j  =  sparsePattern->index[nth];
+            setJacElementKluSparse(j, ii, (newdelta[j] - delta[j]) * delta_hh[ii], nth, Jac);
+            nth++;
+          };
+        }
+        else
+        {
+          nth = (ii == 0) ?  0 : sparsePattern->leadindex[ii-1];
+          while(nth < sparsePattern->leadindex[ii])
+          {
+            j  =  sparsePattern->index[nth];
+            setJacElementKluSparse(j, ii, (newdelta[j] - delta[j]) * delta_hh[ii], nth, Jac);
+            nth++;
+          };
+        }
+        states[ii] = ysave[ii];
+        if (idaData->daeMode)
+        {
+          yprime[ii] = ypsave[ii];
+        }
+      }
     }
   }
   /* finish matrix colptrs */
-  Jac->colptrs[idaData->N] = nth;
-
-  /* not sure if the transposed matrix is needed */
-  /* transposeJac(Jac); */
+  Jac->colptrs[idaData->N] = idaData->NNZ;
 
   unsetContext(data);
 
@@ -1341,7 +1300,7 @@ int jacobianSparseNumIDA(double tt, N_Vector yy, N_Vector yp, N_Vector rr, SlsMa
 }
 
 /*
- * provides a numerical Jacobian to be used with IDA
+ * Wrapper function for jacobianSparseNumIDA
  */
 static int jacobianSparseNum(double tt, double cj,
     N_Vector yy, N_Vector yp, N_Vector rr,
@@ -1367,7 +1326,6 @@ static int jacobianSparseNum(double tt, double cj,
     PrintSparseMat(Jac);
   }
 
-
   /* add cj to diagonal elements and store in Jac */
   if (!idaData->daeMode)
   {
@@ -1385,55 +1343,4 @@ static int jacobianSparseNum(double tt, double cj,
   return 0;
 }
 
-/*
- * provides a numerical Jacobian to be used with IDA
- */
-static int jacobiancoloredKLUNum(double tt, double cj,
-    N_Vector yy, N_Vector yp, N_Vector rr,
-    SlsMat Jac, void *user_data,
-    N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
-{
-  TRACE_PUSH
-
-  IDA_SOLVER* idaData = (IDA_SOLVER*)user_data;
-  DATA* data = (DATA*)(((IDA_USERDATA*)idaData->simData)->data);
-  threadData_t* threadData = (threadData_t*)(((IDA_USERDATA*)((IDA_SOLVER*)user_data)->simData)->threadData);
-
-  SetToZero(idaData->denseJac);
-
-
-  if(jacOwnNumColoredIDA(tt, yy, yp, rr, idaData->denseJac, cj, user_data))
-  {
-    throwStreamPrint(threadData, "Error, can not get Matrix A ");
-    TRACE_POP
-    return 1;
-  }
-  CopySparseMat(SlsConvertDls(idaData->denseJac), Jac);
-
-  /* debug */
-  if (ACTIVE_STREAM(LOG_JAC)){
-    infoStreamPrint(LOG_JAC, 0, "##IDA## Sparse Matrix A.");
-    PrintSparseMat(Jac);
-  }
-
-
-  /* add cj to diagonal elements and store in Jac */
-  if (!idaData->daeMode)
-  {
-    int i;
-    for (i=0; i < idaData->N; ++i){
-      idaData->tmpJac->colptrs[i] = i;
-      idaData->tmpJac->rowvals[i] = i;
-      idaData->tmpJac->data[i] = -cj;
-    }
-    idaData->tmpJac->colptrs[idaData->N] = idaData->N;
-    SlsAddMat(Jac, idaData->tmpJac);
-  }
-
-  TRACE_POP
-  return 0;
-}
-
-
-
 #endif
diff --git a/SimulationRuntime/c/simulation/solver/ida_solver.h b/SimulationRuntime/c/simulation/solver/ida_solver.h
index 1af7bcc0f5b..4810bd8a74b 100644
--- a/SimulationRuntime/c/simulation/solver/ida_solver.h
+++ b/SimulationRuntime/c/simulation/solver/ida_solver.h
@@ -81,6 +81,7 @@ typedef struct IDA_SOLVER
   /* ### daeMode ### */
   int daeMode;                  /* if TRUE then solve dae more with a reals residual function */
   long int N;
+  long int NNZ;
   double *states;
   double *statesDer;
   int (*residualFunction)(double time, N_Vector yy, N_Vector yp, N_Vector res, void* userData);
diff --git a/SimulationRuntime/c/util/simulation_options.c b/SimulationRuntime/c/util/simulation_options.c
index 4eca7149039..4d38c44b5b9 100644
--- a/SimulationRuntime/c/util/simulation_options.c
+++ b/SimulationRuntime/c/util/simulation_options.c
@@ -597,7 +597,6 @@ const char *JACOBIAN_METHOD[JAC_MAX+1] = {
   "numerical",
   "symbolical",
   "kluSparse",
-  "kluColored",
 
   "JAC_MAX"
 };
@@ -611,7 +610,6 @@ const char *JACOBIAN_METHOD_DESC[JAC_MAX+1] = {
   "numerical jacobian.",
   "symbolic jacobian - needs omc compiler flags +generateSymbolicJacobian or +generateSymbolicLinearization.",
   "sparse jacobian for KLU",
-  "colored jacobian for KLU",
 
   "JAC_MAX"
  };
diff --git a/SimulationRuntime/c/util/simulation_options.h b/SimulationRuntime/c/util/simulation_options.h
index 742bedd73a8..6c0417abc78 100644
--- a/SimulationRuntime/c/util/simulation_options.h
+++ b/SimulationRuntime/c/util/simulation_options.h
@@ -241,7 +241,6 @@ enum JACOBIAN_METHOD
   NUMJAC,
   SYMJAC,
   KLUSPARSE,
-  KLUCOLORED,
 
   JAC_MAX
 };