In Some Notes on Internal Implementation especially in Algebra and Calculus one finds interesting subtleties and differences between these two functions, e.g.
The code for Solve and related functions is about 500 pages long.
Reduce and related functions use about 350 pages of Mathematica code
and 1400 pages of C code.
There is much more than a little difference between them. The most important is :
- Reduce returns results of computation as boolean formulae and gives complete description of solution sets.
-
Solve returns lists of replacement rules yielding generic solutions.
An important step toward more complete description of solution sets was a new option of Solve in Mathematica 8, namely MaxExtraConditions (default value 0). Using MaxExtraConditions -> All in Solve provides complete solutions for algebraic equations, nevertheless we have to emphasize that sometimes we might better work with Reduce rather than Solve (regardless of any options added) because replacement rules may appear not a good fit in description of solutions in the real or complex domain to algebraic equations as well as to trancendental equations˛ Distinction between genericity and completness does not make sense in the Integers, an example provided below.
Solve[expr,vars] assumes by default that quantities appearing algebraically in
inequalities are real, while all other quantities are complex.
We have the same issue with Reduce[expr,vars].
Example : Solve cannot find solutions in the real domain
Consider a simple symbolic case in the real domain where Solve does not work even with MaxExtraConditions -> All :
Solve[ x^2 + y^2 <= r^2, {x, y}, MaxExtraConditions -> All]
Solve::fdimc: When parameter values satisfy the condition r ∈ Reals, the solution set
contains a full-dimensional component; use Reduce for complete solution information. >>
{}
Reduce[ x^2 + y^2 <= r^2, {x, y}]
r ∈ Reals && ((x == -Sqrt[r^2] && y == 0) ||
(-Sqrt[r^2] < x < Sqrt[r^2] && -Sqrt[r^2 - x^2] <= y <= Sqrt[r^2 - x^2]) ||
(x == Sqrt[r^2] && y == 0))
We have found none solutions with Solve. We should remember that we can find instances in the real domain with :
fi = FindInstance[ x^2 + y^2 <= 9, {x, y}, 5]
{{x -> 83/84, y -> 31/18}, {x -> 3, y -> 0}, {x -> 83/84, y -> (13 Sqrt[335])/84},
{x -> 5/28, y -> -(Sqrt[7031]/28)}, {x -> 37/21, y -> (10 Sqrt[26])/21}}
However for the same inequality we can find all solutions with Solve in the integers.
Example : Solve provides solutions in the integer domain
s = {x, y} /. Solve[ x^2 + y^2 <= 9, {x, y}, Integers];
r = (List @@@ List @@ Reduce[ x^2 + y^2 <= 9, {x, y}, Integers])[[All, All, 2]]
r == s
True
Here the blue region represents solutions described by Reduce while the red points on the left plot are instances found with FindInstance in the real domain, whereas on the right plot there are all solutions found by Solve in the Integers :
GraphicsRow[{ RegionPlot[ x^2 + y^2 <= 9, {x, -5, 5}, {y, -5, 5},
Epilog -> {Red, PointSize[0.015], Point[fi[[All, All, 2]]]}],
RegionPlot[ x^2 + y^2 <= 9, {x, -5, 5}, {y, -5, 5},
Epilog -> {Red, PointSize[0.015], Point[s]}]}]

Example : Genericity can be advantageous for symbolic transcendental equations
Genericity of Solve output sometimes can be advantageous, e.g. compare
Solve[ x^c == 5, x]
Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be
found; use Reduce for complete solution information. >>
{{x -> 5^(1/c)}}
Solve yields this solution immediately, while Reduce after ~ 10 seconds yields a huge boolean formula, here is a little part of it :
Reduce[ x^c == 5, x][[3, 2]]
C[1] <= -1 && (2 π C[1] - Sqrt[4 π ^2 C[1]^2 + Log[5]^2])/( 2 π ) < Re[c] < 2 C[1] &&
Im[c] <= (-Log[5] - Sqrt[ Log[5]^2 + 8 π ^2 C[1] Re[c] - 4 π ^2 Re[c]^2])/(2 π ) &&
x == E^((2 I π C[1] + Log[5])/c)
There are another examples where we would rather use Reduce rather than Solve with specification MaxExtraConditions -> All because a list of replacement rules cannot express the full solution unlike a boolean form.
Example : Solve cannot find the complete solution set
Even with MaxExtraConditions -> All option Solve sometimes fails which may be seen when we work e.g. with real variables, e.g. :
Normal @ Solve[ a x^2 == c && b x^3 == d, x, Reals, MaxExtraConditions -> All]
Solve::fdimc: When parameter values satisfy the condition d == 0 && b== 0 &&
c == 0 && a == 0,
the solution set contains a full-dimensional component; use Reduce for complete
solution information. >>
{{x -> 0}, {x -> -Sqrt[(c/a)]}, {x -> Sqrt[c/a]}, {x -> Root[-d + b #1^3 &, 1]}}
We have used here Normal to get rid of ConditionalExpressions.
Reduce[ a x^2 == c && b x^3 == d, x, Reals]

For a nice discussion of these topics see this presentation Getting the Most from Algebraic Solvers in Mathematica by Adam Strzeboński during Wolfram Technology Conference 2011.
Example : Further issues
Solve yields exlicit formulae for the solutions (in terms of radicals) to equations whenever it is possible (it is always possible for univariate algebraic (polynomial) equations up to the order four). For more detailed discussion see e.g. : How do I work with Root objects?.
You can find more specific differences examining carefully this :
Options /@ {Reduce, Solve} // Column
{ Backsubstitution -> False, Cubics -> False, GeneratedParameters -> C,
Method -> Automatic, Modulus -> 0, Quartics -> False,
WorkingPrecision -> Infinity}
{ Cubics -> True, GeneratedParameters -> C, InverseFunctions -> Automatic,
MaxExtraConditions -> 0, Method -> Automatic, Modulus -> 0, Quartics -> True,
VerifySolutions -> Automatic, WorkingPrecision -> Infinity }
And even more systematic discussion of the differences might include SystemOptions["ReduceOptions"]. For some sketches of their usage take a look at this reference Real Polynomial Systems.
SystemOptions["ReduceOptions"] // Short[#, 5] &
{ "ReduceOptions" -> {"ADDSolveBound" -> 8, "AlgebraicNumberOutput" -> True,
"BDDEliminate" -> Automatic, "BooleanInstanceMethod" -> Automatic, <<21>>,
"UseNestedRoots" -> Automatic, "UseOldReduce" -> False,
"UseTranscendentalRoots" -> Automatic, "UseTranscendentalSolve" -> True}}
These issues could be systematically discussed only with concrete examples, so that the main question has to be investigated on a case-by-case basis. See e.g. Solving/Reducing equations in $\mathbb{Z}/p\mathbb{Z}$ for a discussion of ExhaustiveSearchMaxPoints.
Edit
Another answer by Itai Seggev says that by specifying Method -> Reduce in Solve, Solve will use Reduce behind the scenes to produce an answer. That is not quite true.
Here is an example : Solve with Method -> Reduce gives a different result than Reduce.
Another statement therein not precise or incorrect is :
"Reduce can deal with the following logical conjunction of inequalities, while Solve cannot."
Solve works well with logical conjunction of inequalities (see e.g. Solve an equation in $\mathbb{R}^+$), let's slightly modify the system :
Solve[ x^2 < 1 && Sin[x] == Cos[x], x]
Reduce[ x^2 < 1 && Sin[x] == Cos[x], x]
Reduce[x^2 < 1 && Sin[x] < Cos[x], x]
{{x -> -2 ArcTan[1 - Sqrt[2]]}}
x == -2 ArcTan[1 - Sqrt[2]]
-1 < x < -2 ArcTan[1 - Sqrt[2]]
whereas the issue here is just mentioned above that the solution set contains a full-dimensional component (in the real domain).
For comparison let's see how it works in the integer domain :
Solve[ x^2 < 2 && Sin[x] < Cos[x], x, Integers]
Reduce[ x^2 < 2 && Sin[x] < Cos[x], x, Integers]
{{x -> -1}, {x -> 0}}
x == -1 || x == 0