Index: ../trunk-jpl/src/c/solvers/solver_newton.cpp
===================================================================
--- ../trunk-jpl/src/c/solvers/solver_newton.cpp	(revision 14166)
+++ ../trunk-jpl/src/c/solvers/solver_newton.cpp	(revision 14167)
@@ -55,12 +55,32 @@
 		xdelete(&old_uf);old_uf=uf;
 
 		/*Solver forward model*/
-		femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL);
+		if(count==1){
+			femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL);
+			CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
+			Reduceloadx(pf,Kfs,ys);xdelete(&Kfs);
+			Solverx(&uf,Kff,pf,old_uf,df,femmodel->parameters);xdelete(&df);
+			Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);xdelete(&ys);
+			InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);xdelete(&ug);
+			xdelete(&old_ug);old_ug=ug;
+			xdelete(&old_uf);old_uf=uf;
+		}
+		uf=old_uf->Duplicate(); old_uf->Copy(uf);
+
+		/*Prepare next iteration using Newton's method*/
+		femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL);xdelete(&df);
 		CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
-		Reduceloadx(pf,Kfs,ys);xdelete(&Kfs);
-		Solverx(&uf,Kff,pf,old_uf,df,femmodel->parameters);xdelete(&df);
+		Reduceloadx(pf,Kfs,ys);   xdelete(&Kfs);
+
+		pJf=pf->Duplicate();
+		Kff->MatMult(uf,pJf);// xdelete(&Kff);
+		pJf->Scale(-1.0); pJf->AXPY(pf,+1.0);     //xdelete(&pf);
+
+		femmodel->CreateJacobianMatrixx(&Jff,kmax);
+		Solverx(&duf,Jff,pJf,NULL,NULL,femmodel->parameters); xdelete(&Jff); xdelete(&pJf);
+		uf->AXPY(duf, 1.0); xdelete(&duf);
 		Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);xdelete(&ys);
-		InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);xdelete(&ug);
+		InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
 
 		/*Check convergence*/
 		convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters); 
@@ -81,20 +101,6 @@
 			break;
 		}
 
-		/*Prepare next iteration using Newton's method*/
-		femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL);
-		CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
-		Reduceloadx(pf,Kfs,ys);   xdelete(&Kfs);
-
-		pJf=pf->Duplicate(); Kff->MatMult(uf,pJf); xdelete(&Kff);
-		pJf->Scale(-1.0); pJf->AXPY(pf,+1.0);     xdelete(&pf);
-
-		femmodel->CreateJacobianMatrixx(&Jff,kmax);
-		Solverx(&duf,Jff,pJf,NULL,NULL,femmodel->parameters); xdelete(&Jff); xdelete(&pJf);
-		uf->AXPY(duf, 1.0); xdelete(&duf);
-		Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);xdelete(&ys);
-		InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
-
 		count++;
 	}