Index: /issm/trunk-jpl/src/c/solvers/solver_newton.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solvers/solver_newton.cpp	(revision 11336)
+++ /issm/trunk-jpl/src/c/solvers/solver_newton.cpp	(revision 11337)
@@ -68,13 +68,15 @@
 
 		/*Prepare next iteration using Newton's method*/
-		SystemMatricesx(&Kff,NULL,&pf,NULL,&kmax,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
-		CreateJacobianMatrixx(&Jff,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,kmax);
-		CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);VecScale(ys,0.);
+		SystemMatricesx(&Kff,&Kfs,&pf,NULL,&kmax,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+		CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
+		Reduceloadx(pf,Kfs,ys);   MatFree(&Kfs);
+
 		VecDuplicate(pf,&pJf);
 		MatMultPatch(Kff,uf,pJf); MatFree(&Kff);
 		VecScale(pJf,-1.);
 		VecAXPY(pJf,+1.,pf);      VecFree(&pf);
-		Solverx(&duf,Jff,pJf,NULL,NULL,femmodel->parameters);
-		MatFree(&Jff);            VecFree(&pJf);
+
+		CreateJacobianMatrixx(&Jff,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,kmax);
+		Solverx(&duf,Jff,pJf,NULL,NULL,femmodel->parameters); MatFree(&Jff);//VecFree(&pJf);
 		VecAXPY(uf,1.,duf);      VecFree(&duf);
 		Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);VecFree(&ys);
Index: /issm/trunk-jpl/src/m/solvers/solver_newton.m
===================================================================
--- /issm/trunk-jpl/src/m/solvers/solver_newton.m	(revision 11336)
+++ /issm/trunk-jpl/src/m/solvers/solver_newton.m	(revision 11337)
@@ -51,5 +51,4 @@
 		Jff=CreateJacobianMatrix(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,kmax);
 		duf=Solver(Jff,pJf,[],[],femmodel.parameters);
-		ys=0*ys;
 		uf=uf+duf;
 		ug=Mergesolutionfromftog(uf,ys,femmodel.nodes,femmodel.parameters);