What happens to Newtonian systems as the mass vanishes? This question is closely related to another one I asked recently, and may be thought of as a warm-up to that one.
Consider $\mathbb R^n$ with its usual metric, and pick a one-form $b$ and a function $c$.  Let $m$ be a positive constant, and consider the second-order differential equation for a function $q(t)$
$$ m\ddot q = db \cdot \dot q + dc $$
where I have used the metric to identify vectors and covectors, $dc$ is the differential of $c$, and $db$ is the exterior derivative of $b$ (it is contracted with $\dot q$ to yield a covector).  In coordinates, and using Einstein's summation convention:
$$ m\ddot q^i = \left(\partial\_i b\_j - \partial\_j b\_i\right)\dot q^j + \partial\_i c $$
I am interested in the limit when $m\to 0$.  For example, when $m=0$ and $b=0$ (or anyway when $b$ is closed), then the differential equation forces the path $q(t)$ to stay within the set of critical points of $c$ (this set is generically discrete, so that the only solutions are constant).  At another (more generic) extreme, $db$ might be nondegenerate, and hence a symplectic form on $\mathbb R^n$.  Then the equation $0 = db \cdot \dot q + dc$ is a nondegenerate first-order differential equation, exactly equivalent to Hamilton's equations for the symplectic manifold $(\mathbb R^n,db)$ with Hamiltonian $-c$.  There is some gradation when $db$ is nonzero but has nontrivial kernel (as for example must happen if $n$ is odd).
So I basically get what happens when $m=0$.  But can we understand the limit $m\to 0$?  For example, if $m\neq 0$, then any initial value $(\dot q(0),q(0))$ determines a solution; for fixed initial values, how does this solution vary as $m\to 0$?  Alternately, we can try to solve the boundary value problem, in which we prescribe $q(0)$ and $q(1)$.  Then what happens to the solutions as $m$ shrinks?  Since when $m=0$ we cannot find solutions with arbitrary initial velocity, it is unlikely that anything is particularly well-behaved in the limit, but not impossible.
Very specifically, I would like to know about the asymptotics of the solutions to the boundary and initial-value problems — what do solutions look like when $m$ is a formal variable?  But more generally I'm happy with some statements about the regularity in the $m\to 0$ limit.
 A: This is not really an answer to your question, but a philosophical point.  The physics of massless particles is not simply the limit $m\to 0$ of a massive particle.  This is perhaps easier to see in the context of relativistic field theory, where free field equations can be interpreted in terms of the representation theory of the Poincaré group.  Let's consider a massive "vector" particle in 4-dimensional Minkowski spacetime.  This translates into the irreducible unitary representation of the Poincaré group which is induced (à la Wigner, Mackey,...) from the three-dimensional representation of the stability subgroup of a momentum with nonzero mass, which is isomorphic to SO(3) for all nonzero values of the mass.  The fact that the inducing representation is three-dimensional explains the physical statement that massive vectors have three degrees of freedom.  
Massless vectors, on the other hand, have only two physical degrees of freedom: the two transversal polarisations of light.  The reason is that they are induced from a real two-dimensional representation of the maximal compact subgroup of the stability subgroup of a nonzero momentum with zero mass, which is isomorphic to SO(2).
The massive and massless cases are thus very different and you cannot view one as a limit of the other, at least from the point of view of the representation theory.  Physically, what is going on is that in the limit $m\to 0$, one of the physical polarisations of the vector becomes 'gauge' and hence unphysical.
In the lagrangian formulation you can see this very clearly, since the mass enters multiplying the Minkowski norm of the vector field $\tfrac12 m \int |A|^2$ and you can easily set it to zero, but then you see that the lagrangian becomes degenerate, signaling that you have a constrained system,...
A: These are called singularly perturbed problems, and in general very large gradients called boundary layers will develop in the solution. For instance, in the boundary value problem you mentioned, as mass goes to zero, the boundary condition at t=1 becomes more and more "invisible" to the bulk of the solution, so that it will closely resemble the solution of the initial value problem with mass=0, until it reaches very close to t=1, where it rushes towards q(1)=0 to satisfy the boundary condition. For the initial value problem with very small mass, similar layer will form near t=0.
A: This is not a real answer, but there is a branch of mathematics called "semiclassical analysis" which might be related. For example, consider a degenerate version of the problem above:
$$
(-h^2 \partial_x^2+V(x))u=0.
$$
Here $h^2=m$ and $V=dc$; we assume that $n=1$. Then the limit as $h\to 0$ is called "semiclassical limit". What should happen (if you prescribe some boundary conditions) is that the possible solutions $u$ should get "microlocalized" near the zero set $\{p=0\}$ of the semiclassical symbol
$$
p(x,\xi)=\xi^2+V(x).
$$
Here a function $u$ is "microlocalized" near a subset $K$ of the cotangent bundle if a certain norm $\|Au\|$ is small for any pseudodifferential operator $A$ with symbol $a$ supported outside of $K$.
A physical explanation of the above would be that our static Schr\"odinger equation governs the behaviour of a single quantum particle under the potential $V$ near the zero energy level; for small values of Planck constant $h$, this should correspond to the motion of a classical particle at this fixed energy level.
There are several sources to read about semiclassical analysis, including 
lecture notes by Evans-Zworski and a book by Dimassi and Sjostrand.
