Force between two point dipoles What would the force between two parallel point dipoles be?
I was thinking of doing it the way force between two point charges is found out,   by finding the field and then the force but I am not able to formulate it.
 A: Going from the field to the force will be difficult, because you don't know a priori how a dipole reacts to a field, particularly if the field is not homogeneous. The way to do this is the same as you find the field for the first dipole: find the net force on the two charges of the second dipole, and then take the limit as their separation goes to zero (whilst keeping the charge large enough so the dipole moment is constant). This is, in fact, rather messy, so it is better to compute the interaction energy, and then find the force as the gradient of this energy.
To do that, start with a point dipole with dipole moment $\mathbf p$ at the origin, which causes an electrostatic potential of
$$V_\text{dip}=\frac{\mathbf p\cdot \mathbf r}{r^3},$$
so the potential energy of a charge $q$ at position $\mathbf r$ is 
$$U_{q}=q\frac{\mathbf p\cdot \mathbf r}{r^3}.$$
To make the second dipole, start with a charge $-q$ at position $\mathbf r$, and add a second charge $+q$ at position $\mathbf r+\Delta r\,\hat{\mathbf{n}}$. The dipole moment of the pair is $\mathbf p'=q\Delta r\,\hat{\mathbf{n}}$, and their interaction energy with the first dipole is
$$U_\text{fin. dip.}=-q\frac{\mathbf p\cdot \mathbf r}{r^3}+q\frac{\mathbf p\cdot(\mathbf r+\Delta r\,\hat{\mathbf{n}})}{\|\mathbf r+\Delta r\,\hat{\mathbf{n}}\|^3}.$$
Now comes the tricky business of taking the limit $\Delta r\to 0$. What makes it complicated is the presence of $\Delta r$ in the denominator, which must be dealt with via the first two terms of a binomial series:
$$\begin{align}
\|\mathbf r+\Delta r\,\hat{\mathbf{n}}\|^{-3}&=\left(r^2+2\Delta r\,\hat{\mathbf{n}}\cdot\mathbf r +\Delta r^2\right)^{-3/2}
\\& =\frac{1}{r^3}\left(1-\frac32 \times \frac{2\Delta r\,\hat{\mathbf{n}}\cdot\mathbf r}{r^2}+O\left(\frac{\Delta r^2}{r^2}\right)\right)
\\& \approx\frac{1}{r^3}-3\frac{\Delta r}{r^4}\hat{\mathbf{n}}\cdot\hat{\mathbf{r}},
\end{align}$$
ignoring the quadratic terms. Putting this into the interaction energy, you get
$$U_\text{fin. dip.}=-q\frac{\mathbf p\cdot \mathbf r}{r^3}+q(\mathbf p\cdot\mathbf r+\Delta r\,\mathbf p\cdot\hat{\mathbf{n}})\left(\frac{1}{r^3}-3\frac{\Delta r}{r^4}\hat{\mathbf{n}}\cdot\hat{\mathbf{r}} \right)+O\left(q\frac{\Delta r^2}{r^5}\right).$$
Here the constant terms cancel, and you're left with two different linear terms.
$$U_\text{fin. dip.}=
\frac{(q\Delta r\,\hat{\mathbf{n}})\cdot\mathbf p}{r^3}
-3\frac{1}{r^3}(q\Delta r\hat{\mathbf{n}})\cdot\hat{\mathbf{r}}\times\mathbf p\cdot\hat{\mathbf{r}}
+O\left(q\frac{\Delta r^2}{r^5}\right).$$
To take the limit, you can note that in the linear terms stay constant at the dipole moment $\mathbf p'=q\Delta r\,\hat{\mathbf{n}}$, so they just stay as they are. The quadratic terms, on the other hand, have a charge that increases as $q=p'/\Delta r$, but their $\Delta r$ dependence is faster and they drop to zero. In the limit, then, the interaction energy is
$$U_\text{dip.-dip.}=
\frac{
\mathbf p'\cdot\mathbf p
-3(\mathbf p'\cdot\hat{\mathbf{r}})(\mathbf p\cdot\hat{\mathbf{r}})
}{r^3}
.$$
To get the force felt on the second dipole, you need to take the gradient with respect to $\mathbf r$, or repeat this calculation with the dipole force
$$
\mathbf F_\text{dip}=-\nabla V_\text{dip}=-\frac{\mathbf p}{r^3}+3\frac{\mathbf p \cdot \mathbf r}{r^5} \mathbf r.
$$
The former is simpler, though you may want to use the form 
$$U_\text{dip.-dip.}=
\frac{
\mathbf p'\cdot\mathbf p
}{r^3}
-\frac{
3(\mathbf p'\cdot{\mathbf{r}})(\mathbf p\cdot{\mathbf{r}})
}{r^5}
.$$
Nevertheless, unless you can fully reproduce the steps in my calculation, I'd recommend you repeat it with the force, and check that both approaches coincide, or there will be something you're missing.
