void SpringNetwork::Simulate(float3N &Xnew, float3N &Vnew, float dt, ClothVertex *extforce)
{	
	// Get ready for conjugate gradient iterative solver step.
	// Initialize operands.
		
	int n=X.count;  // all our big vectors are of this size
	float3N temp(n);  // temp to store result of matrix multiply
	float3N B(n);  
	
	dV.Zero();	
	A.Identity();  // build up the big matrix we feed to solver
		
	if(!step) // implicit euler as initial guess for 2-nd order bdf
	{	
		A = M * A;
		A -= (dFdV * dt +  dFdX * (dt*dt)); // A = I - (dFdV * dt +  dFdX * (dt*dt))
		
		temp = dFdX * V;
		B = F * dt + temp * (dt*dt);
		
	}
	else // BDF-2
	{                                           		 
		
		// multiplying the whole equation from left with M
		// idea taken from Baraff[98] who applied this to the implicit euler formula
		
		A = M * A; // M - [...]
		A -= (dFdV * (2.0f*dt/3.0f) + dFdX * ((4.0f*dt*dt)/9.0f)); // A = I - (dFdV * (2*dt/3) +  dFdX * (4*dt*dt/9))
				
		B = (X-Xold) * (1.0f/3.0f); // (X - Xold) * 1/3 	
		
		B += ( (V*8.0f -  Vold*2.0f) * (dt/9.0f)); // + (8*V - 2*Vold)*(dt/9);
		
		B = M * B; 
		
		temp = dFdV * V;  // + ((4*dt*dt)/9)*(F - dFdV*V)
		B += (F - temp) * ((4.0f*dt*dt)/9.0f); 
		
		temp = (X - Xold);  // - (2*dt/9)*dFdV*(X-Xold);
		temp = dFdV*temp;
		B = B - (temp * (2.0f*dt/9.0f)); 	
	}	
	
	PreConjGradientFiltered(dV,A,B,S,z, olddV);	
	
	olddV = dV;
	
	if(!step)
	{	
		Vnew = V + dV;
		Xnew = X + V*dt; 
		
	}
	else
	{
		Xnew = X + dV;
		
		Vnew = Xnew*(3.0f/2.0f) - X*2.0f;
		Vnew = (Vnew + Xold*0.5f);
		Vnew = Vnew * dt;
	}
}