Index: /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp	(revision 11323)
+++ /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp	(revision 11324)
@@ -143,7 +143,11 @@
 	/*Assign output pointers: */
 	if(pKff) *pKff=Kff;
+	else      MatFree(&Kff);
 	if(pKfs) *pKfs=Kfs;
+	else      MatFree(&Kfs);
 	if(ppf)  *ppf=pf;
+	else      VecFree(&pf);
 	if(pdf)  *pdf=df;
+	else      VecFree(&df);
 	if(pkmax) *pkmax=kmax;
 }
Index: /issm/trunk-jpl/src/c/solvers/solver_newton.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solvers/solver_newton.cpp	(revision 11323)
+++ /issm/trunk-jpl/src/c/solvers/solver_newton.cpp	(revision 11324)
@@ -47,31 +47,18 @@
 	for(;;){
 
-		//save pointer to old velocity
 		VecFree(&old_ug);old_ug=ug;
 		VecFree(&old_uf);old_uf=uf;
 
+		/*Solver forward model*/
 		SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
 		CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
 		Reduceloadx(pf,Kfs,ys);MatFree(&Kfs);
-		Solverx(&uf,Kff,pf,old_uf,df,femmodel->parameters);
+		Solverx(&uf,Kff,pf,old_uf,df,femmodel->parameters);VecFree(&df);
+		Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);VecFree(&ys);
+		InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);VecFree(&ug);
 
-		convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters);
-
-		CreateJacobianMatrixx(&Jff,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
-
-		Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);
-		InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
-		SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
-		VecDuplicate(pf,&pJf);
-		MatMultPatch(Kff,uf,pJf);
-		VecAXPY(pJf,-1.,pf);
-		Solverx(&duf,Jff,pJf,NULL,NULL,femmodel->parameters);
-		MatFree(&Kff);VecFree(&pf); VecFree(&df);
-		VecAXPY(uf,-1.,duf);
-		Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);VecFree(&ys);
-
-		InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
-
-		/*Increase count: */
+		/*Check convergence*/
+		convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters); 
+		MatFree(&Kff);VecFree(&pf);
 		count++;
 		if(converged==true)break;
@@ -80,4 +67,18 @@
 			break;
 		}
+
+		/*Prepare next iteration using Newton's method*/
+		CreateJacobianMatrixx(&Jff,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+		SystemMatricesx(&Kff,NULL,&pf,NULL,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+		CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
+		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);
+		VecAXPY(uf,+1.,duf);      VecFree(&duf);
+		Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);VecFree(&ys);
+		InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
 	}