First, we will need to determine a form for the function $f(r)$ which describes the brightness at a distance $r$ from the center of the star. This is the effect of limb darkening; there are several models. We can implement a very rudimentary one geometrically by finding the angle between a surface normal and a ray cast toward us from that point on the surface (at a projected distance of $r$). See the following:
From this figure, you can see that for a surface normal the component directed toward the observer will be $\cos\theta$ (taking $1$ as the ``length'' of the initial ray), where $\theta$ is the angle between the surface normal and the observer. Likewise, we can see that the angle depends on the projected distance $r$ such that $\theta=\sin^{-1} r$. So, we can finally write our brightness function as
$$f(r)=\begin{cases}
A\cos\sin^{-1}r&r\leq1\\0&r\geq1
\end{cases}$$where $A$ is a normalization constant. We will want to use it to ensure that the total brightness of the unocculted star is $1$. To that end, we can determine $A$ as such: $$\begin{align*}
1&=\int_0^1\int_0^{2\pi}\,\mathrm{d}\theta\,\mathrm{d}r\,Ar\cos\sin^{-1}r\\
1&=\frac{2\pi A}{3}\\
\implies A&=\frac{3}{2\pi}
\end{align*}$$
This brightness profile looks like follows:
Consider that we can uniquely define the track across the star using two parameters, the radius of the secondary body $R$, and the impact parameter $b$. Here, $b$ just gives the projection distance closest approach between the centers of the two bodies. If we choose coordinates such that the track of the transiting body is horizontal (along the $x$-coordinate), then we can describe the line of the track as $y=b$. This choice of coordinates also allows us to parameterize the motion completely along the $x$ direction.
The brightness $B(t)$ at a position along the track $x=t$ (coordinate of the center of the occulting body) for some $b$ and $R$ can then be determined by integrating the brightness function over the area of the disk of the secondary, and subtracting it from $1$:
$$\begin{align*}
B(t)&=1-\int_D\mathrm{d}A\,f(r)\\
&=1-\frac{3}{2\pi}\int_{t-R}^{t+R}\int_{b-\sqrt{R^2-(x-t)^2}}^{b+\sqrt{R^2-(x-t)^2}}\mathrm{d}y\,\mathrm{d}x\,f\left(\sqrt{x^2+y^2}\right)
\end{align*}$$
This form is for an occulting body which is spherical, totally opaque, has no surrounding dust, does not change in apparent diameter over the duration of the transit, etc. It is general in terms of the brightness profile - if you want to use a better model (for which I'd recommend doing some reading, perhaps starting here) you can simply replace $f(r)$. Anyway, I don't think it is likely one will be able to find some analytical form for this in general, without making some simplifying assumptions (for example, we could assume that the occulting body is sufficiently small that the brightness of the star does not vary significantly over its area). It lends itself nicely to numerical integration, though.
Here is an animation showing 20 stages along a transit using this method. The lightcurve plot is generated from 100 samples. This particular situation is for $b=0.1$, $R=0.02$.
Compare against some Kepler data (source):
