How can I limit the amount of curvature of a bezier curve?

First of all, you need to match the endpoints, as you said. I will assume that you also need to match the tangents (especially if it's a part of a spline). If the control points are A,B,C,D, the tangent in A points along vector AB, and the tangent in D points in DC direction.

You can therefore freely move B and C along the tangent direction, which changes the "sharpness" but keeps the tangents and endpoints intact.

There are three typical cases:

1) The tangents converge. This means that intersection of lines AB and CD is where the curve is. The closer B and C are to this intersection, the sharper the curve. If B or C is farther away than the intersection, you have a self-intersecting curve or at least something extremely sharp. If your condition is "soft" (not mathematical, but you just want some heuristic to limit the curvature), you can for instance, just forbid B and C to be more than half way from A/D to the intersection. If you want mathematically rigorous condition, you can vary the two parameters (AB and CD distances) and calculate the minimal radius of curvature for each case until you go over the desired limit.

2) the tangents diverge on the same side. This case produces a "bubble". The farther away the B and C are, the larger the circle-like bump is. This kind of a spline is usually not something you want (at least in graphics), you have more control if you split it in half and get two "conventional" beziers. In any case, here, the curve is sharprer if AB and CD are smaller, but varying these parameters makes drastic changes in the shape.

3) the tangents point in opposite directions (the curve intersects the A-D straight line). This is a horrible case (S-shape). Still, you can only vary AB and CD distances if you want to keep the tangents, so you can vary them until you honor the condition. This case again has subcases (like the two above), depending on whether the angles BAD and ADC are obtuse or acute.

I hope this helps. Essentially, you have two free parameters that you can vary in order to minimize sharpness. It also makes sense (to preserve the general shape) to keep the ratio AB/CD fixed and scale the distances simultaneously. Then you have a 1-parametric case which you can solve algorithmically without ambiguity (in 2-parametric case, you have too many solutions).

---

Edit: I was hoping to avoid explicit math because it's easy to know what to do, but it's tedious to write down and do anything analytically, better to let the computer do it. However, here we go:

A cubic bezier is parameterized like this: $$\vec{r}(t)=(1-t)^3 A + 3(1-t)^2 t B+ 3(1-t)t^2 C +t^3 D$$ the first derivative: $$\vec{r}'(t)=3(1-t)^2(B-A)+6(1-t)t(C-B)+3t^2(D-C)$$ the second derivative is $$\vec{r}''(t)=6(1-t)(C-2B+A)+6t(D-2C+B)$$

Now, the curvature is expressed as such: $$\kappa=\frac{|\vec{r}'\times\vec{r}''|}{|\vec{r}'|^{3}}$$

You pretty much need to calculate this numerically and also find the maximum over $t$ numerically. It's too ugly to do anything "on paper". See the discussion here:

maximum curvature of 2D Cubic Bezier

My suggestion was to "correct" your control points as such: $$B_u=A+(B-A)u$$ $$C_u=D+(C-D)u$$ where you are looking for $u$ which minimizes your maximum curvature. $u=1$ is your original curve. This parametrization lets you keep the approximate shape of the curve.

If you would like to minimize the "sharpness" of the cubic Bezier curve honoring two end points and two end tangent directions, there is something called "optimized geometric Hermite curve" (OGH curve) that might be interested to you. The OGH curve does not minimize the maximum curvature of the curve. Instead, it minimizes the overall "strain energy" of the curve, which is

$$\int_0^1 [f^"(t)]^2 \;\mathrm{dt}$$

You can refer to this paper link for details. For cubic OGH curve, you can find out the "optimal" magnitudes of the end derivatives analytically. The formula is listed in this paper as equation (4).