Index: /issm/trunk-jpl/src/c/solvers/solver_newton.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solvers/solver_newton.cpp	(revision 14415)
+++ /issm/trunk-jpl/src/c/solvers/solver_newton.cpp	(revision 14416)
@@ -15,5 +15,5 @@
 	/*intermediary: */
 	bool   converged;
-	int    count;
+	int    count,newton;
 	IssmDouble kmax;
 	Matrix<IssmDouble>* Kff = NULL;
@@ -37,4 +37,5 @@
 	femmodel->parameters->FindParam(&max_nonlinear_iterations,DiagnosticMaxiterEnum);
 	femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+	femmodel->parameters->FindParam(&newton,DiagnosticIsnewtonEnum);
 	femmodel->UpdateConstraintsx();
 
@@ -56,10 +57,30 @@
 
 		/*Solver forward model*/
-		femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL);
+		if(count==1 || newton==2){
+			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);
+
+		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*/
@@ -82,18 +103,4 @@
 		}
 
-		/*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++;
 	}