function [u_g varargout]=diagnostic_core_nonlinear(m,inputs)
%INPUT function [ru_g varargout]=cielodiagnostic_core_nonlinear(m,inputs)
	
	%check inputs
	if isempty(inputs),
		clear inputs;
	end

%   first off! We are going to modify the loads dataset. We need to shield the loads from those changes once we return;
	loads=m.loads;

%	stiffness and load generation only:
	m.parameters.kflag=1; m.parameters.pflag=1;

%   initialize solution vector
	count=1;
	converged=0;
	soln(count).u_g=[];
	soln(count).u_f=[];


	if m.parameters.debug,
		disp(sprintf('\n%s',['   starting direct shooting method']));
	end
	
	while(~converged),
		
		%add velocity to inputs
		inputs.velocity=full(soln(count).u_g);

		%Update inputs in datasets
		[m.elements,m.nodes, loads,m.materials]=UpdateFromInputs(m.elements,m.nodes, loads,m.materials,inputs);

		%system matrices 
		[K_gg_nopenalty , p_g_nopenalty]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs);
		
		%penalties
		[K_gg , p_g]=PenaltySystemMatrices(K_gg_nopenalty,p_g_nopenalty,m.elements,m.nodes,loads,m.materials,m.parameters,inputs);


		%Reduce tangent matrix from g size to f size
		[K_ff, K_fs] = Reducematrixfromgtof( K_gg, m.Gmn, m.nodesets); 

		%Reduce load from g size to f size
		[p_f] = Reduceloadfromgtof( p_g, m.Gmn, K_fs, m.ys, m.nodesets);

		%Increment index
		count=count+1;

		%Solve	
		[soln(count).u_f]=Solver(K_ff,p_f,[],m.parameters);
	
		%Merge back to g set
		[soln(count).u_g]= Mergesolutionfromftog( soln(count).u_f, m.Gmn, m.ys, m.nodesets ); 

		if (count>2),
			soln(count-1).u_f=NaN;
		end

		%Deal with penalty loads
		inputs.velocity=soln(count).u_g;
	
		%penalty constraints
		[loads,constraints_converged,num_unstable_constraints] =PenaltyConstraints( m.elements,m.nodes, loads, m.materials,m.parameters,inputs);

	%   Figure out if convergence is reached.
		if(count>=3),
			dug=soln(count).u_g-soln(count-1).u_g; ndu=norm(dug,2); nu=norm(soln(count-1).u_g,2); 
			if (ndu/nu<=m.parameters.eps_rel),
				if constraints_converged,
					converged=1;
				end
			end
			
			if m.parameters.debug,
				disp(sprintf('%s %g %s %g','      convergence criterion: norm(du)/norm(u)=',ndu/nu,' > ',m.parameters.eps_rel));
				disp(sprintf('%s %i','      number of unstable constraints=',num_unstable_constraints));
			end
		end

		if (count>2),
			soln(count-1).u_g=NaN;
		end
	end
			
	%some cleanup
	[soln(count).u_f]=NaN;

	u_g=soln(count).u_g;

	%more output might be needed, when running in cielocontrol.m
	nout=max(nargout,1)-1;
	if nout==2,
		
		inputs.velocity=full(soln(count).u_g);
		m.parameters.kflag=1; m.parameters.pflag=0; 
		[K_gg, p_g]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs);
		[K_ff, K_fs] = Reducematrixfromgtof( K_gg, m.Gmn, m.nodesets); 
		varargout(1)={K_ff};
		varargout(2)={K_fs};
	end
	
end
