function femmodel=solver_newton(femmodel)
%SOLVER_NEWTON - core solver of diagnostic run
%
%   Usage:
%      [femmodel]=solver_newton(femmodel)

	%Branch on partitioning schema requested
	maxiter=femmodel.parameters.DiagnosticMaxiter;
	configuration_type=femmodel.parameters.ConfigurationType;
	[femmodel.nodes]=UpdateConstraints(femmodel.nodes,femmodel.constraints,femmodel.parameters);

	%initialize solution vector
	converged=0; count=1;

	%Start non-linear iteration using input velocity: 
	ug=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, femmodel.loads, femmodel.materials, femmodel.parameters);
	uf=Reducevectorgtof( ug, femmodel.nodes,femmodel.parameters);

	%Update the solution to make sure that vx and vxold are similar
	[femmodel.elements femmodel.loads]=InputUpdateFromConstant(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,double(converged),ConvergedEnum);
	[femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ug);

	while(~converged),

		%save pointer to old velocity
		old_ug=ug;
		old_uf=uf;

		%Solver forward model
		[K_ff,K_fs,p_f,df,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
		ys=CreateNodalConstraints(femmodel.nodes,configuration_type);
		p_f=Reduceload( p_f, K_fs, ys);
		issmprintf(VerboseSolver(),'%s%g','      condition number of stiffness matrix: ',condest(K_ff));
		uf=Solver(K_ff,p_f,old_uf,df,femmodel.parameters);
		ug=Mergesolutionfromftog( uf, ys, femmodel.nodes,femmodel.parameters); 
		[femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ug);

		%Check convergence
		converged=convergence(K_ff,p_f,uf,old_uf,femmodel.parameters);
		if(converged==1) break; end
		if(count>maxiter),
			issmprintf(true,'%s%i%s','      maximum number of iterations ',maxiter,' exceeded');
			break;
		end

		%Prepare next iteration using Newton's method
		[K_ff,K_fs,p_f,df,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
		ys=CreateNodalConstraints(femmodel.nodes,configuration_type);
		p_f=Reduceload( p_f, K_fs, ys);
		pJf=p_f-K_ff*uf;
		Jff=CreateJacobianMatrix(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,kmax);
		duf=Solver(Jff,pJf,[],[],femmodel.parameters);
		uf=uf+duf;
		ug=Mergesolutionfromftog(uf,ys,femmodel.nodes,femmodel.parameters); 
		[femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ug);

		%increase count
		count=count+1;
	end
end
