Searching for three positive integers as a solution to an equation
Here are some simple optimizations you can make that can push your limit a bit higher:


Don't use floating point arithmetic. Floating point arithmetic is slow compared to integer arithmetic (this is not always true with the advent of things like FPUs, but it is a good rule of thumb). Even worse, floating point arithmetic is prone to things like rounding errors, which means you might find solutions which don't work and miss solutions that do work.

One way to avoid using floating point arithmetic here is to expand out the equation into the form:

\$x(x+y)(x+z) + y(y+x)(y+z) + z(z+x)(z+y) = 4(y+z)(x+z)(x+y)\$

In this form you can check whether the equation is satisfied just using integer arithmetic (make sure you use a integer type large enough so that neither side overflows, however).
Make use of symmetry (thanks @Deduplicator). Another simple observation you can make is that this equation is symmetric in \$x\$, \$y\$, and \$z\$, i.e. if you have a solution \$(x, y, z)\$, then any permutation of this solution also works. One way to use this fact is to only loop over \$1 \leq y \leq x\$ and \$1 \leq z \leq y\$. This cuts down on your total number of iterations by a factor of \$3! = 6\$. 
Reduce number of degrees of freedom. One interesting observation is that, once you fix \$x\$ and \$y\$, you don't have many choices left for \$z\$. So instead of looping over all \$z\$ in the range \$[1, L]\$, you can more efficiently compute the possible values of \$z\$. (This removes one of your for loops).

How do you find the possible values of \$z\$? Well, if we know what \$x\$ and \$y\$ are, the above equation reduces to a cubic in \$z\$; i.e., something for the form \$c_3z^3 + c_2z^2 + c_1z + c_0 = 0\$. There are a couple ways to solve this for \$z\$; the easiest way I can think of (but definitely not the fastest) is to try all the factors of \$c_0\$ (by the Rational Root Theorem, any integer \$z\$ that satisfies this must divide \$c_0\$). This does involve factoring \$c_0\$, however, which can take a while depending on how you implement it. 

A faster method is to use binary search to find the roots, but you have to be a bit careful here, since the cubic might not be monotone increasing/decreasing over the entire interval. The correct way to do this is to first find the roots of the derivative of the cubic (the derivative of the cubic is a quadratic, so you can do this with just the quadratic formula); these are the critical points of the cubic, so the cubic is monotone increasing/decreasing on the intervals between these points, and you can binary search here. 
Learn about elliptic curves. So it turns out that the smallest solutions to this Diophantine equation are immense, and it's pretty much impossible to find any via any method like this (I forget how large they are, but I believe each of x, y, and z is at least 10 digits long). The question then becomes, how did people find this answer in the first place?

The trick here is that the cubic above (in fact, pretty much any homogeneous cubic in three variables) is an instance of a mathematical object known as an elliptic curve, and that these objects have a ton of nice properties. One of these nice properties is that there is a method whereby you can take any two integer solutions to this equation and "add" them to get a third integer solution. So one approach which works is to start with an integer solution (but not a positive integer solution) such as (11, 9, -5), and repeatedly "add" it to itself until you end up with a positive integer solution (as far as I know, there's no guarantee you ever will end up with a positive integer solution, but in this case you do after a couple of steps). 


EDIT:

I ran some tests to compare the above optimizations (see here). On my machine for limit=1000:


OP's code (removing the printf per iteration) takes ~14.5 seconds.
OP's code with optimization 1 takes ~4.2 seconds.
OP's code with optimizations 1 + 2 takes ~0.7 seconds.
OP's code with optimizations 1 + 2 + 3 takes ~0.25 seconds. 


For larger limits, these differences are more pronounced. For limit=5000:


OP's code with optimizations 1 + 2 + 3 takes ~7.9 seconds.
OP's code with optimizations 1 + 2 takes ~86 seconds.
The other two cases each took at least 5 minutes.