How to integrate over a timestep in a mass-spring model? I'm writing a simulation of a block of matter using a "mass-spring" model, where the matter is modelled as a 3D lattice of point masses, where each point is connected by springs to the (up to) 26 neighboring point masses.
At each "tick" I calculate the new force at each point mass with (pseudocode):
for each neighbor:
    distance = length(vector_from_me_to_neighbor)
    force_here += spring_constant * (distance - ideal_distance_from_me_to_neighbor) * normalised_vector_from_me_to_neighbor

Using verlet integration, I update the position of each point mass with:
new_pos = current_pos + velocity_approximation + (force / point_mass) * timestep * timestep * 0.5;

(The timestep is between 16 and 32ms.)
This works reasonably for very "squishy" materials (with a low spring constant), but for stiffer materials the simulation explodes as the course timestep calculation bounces the masses back and forth further and further.
My question is, is there a resonably efficient way of integrating over the timestep rather than calculating the force once and applying it for the duration of the timestep?  Could I potentially use the first two terms of a taylor expansion?
 A: It sounds like you're running into the problem of ODE stiffness, you have a system in which vastly different timescales need to be considered. There is an entire field of applied mathematics dedicated to the study of numerical ODE solvers, and it's not possible to discern exactly which method would suit you best from the information you've provided.
Your solver is what's known as an explicit method - the $n$th time step is calculated based solely on previous timesteps. Though simple to implement, such methods are intrinsically unstable when used with large time steps. It sounds like you need an implicit method instead. Look into something like the backwards Euler method. Note that this is not symplectic, so it does not conserve energy exactly, nor is it particularly accurate, but it is extremely stable.
A: I suggest that before writing pseudocode, you try to master the underlying algorithm. If what you spelled verlot is supposed to be the Verlet algorithm, you miss the key ingredient that is the approximation for velocities. Moreover, the time-step should appear in your equation for new_pos.
Once you can implement Verlet's algorithm in the velocity form, you may go quite far before you require implicit methods. In any case, do not waste your time with the backward Euler method. It is first-order, while Verlet is second order. Moreover, Verlet's algorithm is a symplectic method. In the numerical simulation of hamiltonian systems, this is a valuable bonus.
