Ok, I'm back with results!
I tried two approaches:
Using solids mechanics to derive a differential equation governing the movement of the centers of the wheels: the inputs of the system "bike" are the torque at the rear wheel and the angle of the front wheel, and the outputs are the kinematics of the centers of the wheels. But I gave up, it was hard!
Trying to guess what happens from a geometric point of view when the rear wheel "pushes" the front wheel forward with the front wheel not straight. This method directly yields an equation of infinitesimal increments (see below) from which you can get an actual differential equation. I haven't tried manipulating this first equation to get the ODE but my guess is that I would have obtained that very same ODE using solids mechanics. It just feels right.
Notations and hypotheses:
We are in the plane with basis vectors ex and ey.
A is the center of the rear wheel. B is the center of the front wheel. The length of the bike L is the distance between A and B. The angle between ey and the vector AB is φ. The angle between AB and the front wheel is θ.
We suppose that, at a certain instant t, A(t) has a velocity V(t) colinear with AB. Therefore, for an infinitesimal timestep dt,
A(t+dt) = A(t) + V(t).dt.
We also suppose that, at time t, the front wheel doesn't drift, i.e. the speed of B is colinear with the direction of the front wheel, i.e. forms an angle θ with AB. We call Uθ the unit vector forming an angle θ with AB, i.e. the unit vector with the same direction as the front wheel.
Therefore, at t+dt,
B(t+dt) = B(t) + λ.Uθ
for a certain real, positive λ such that the length of the bike L is conserved:
distance( A(t+dt) , B(t+dt) ) = L
This last equation translates to
norm²( B(t) + λ.Uθ - A(t) - V(t).dt ) = L²
but B(t), by definition, is A(t) + L.Uφ, so that λ must satisfy the equation
norm²( L.Uφ + λ.Uθ - V(t).dt ) = L².
The solution, of course, is independent from φ since the problem is the same when the bike points towards positive y. Therefore, if we call R the rotation matrix with angle -φ, λ must be the positive solution of
norm²( L.ey; + λ.Uθ - R.V(t).dt ) = L².
After a few calculations, if we call v the norm of V, you get
λ = L.( sqrt( 1 - (sin(θ).(1-v.dt/L))² ) - cos(θ) ) + v.dt.cos(θ).
Here's the pseudocode I used to get the animation above (instead of using Uθ, I use u = U(θ+φ) because it was simpler):
// I start at i=1 because i=0 contains the initial values
for (int i=1; i<=N; i++)
// the array in which I stored the successive A points
Aarray[i] = Aarray[i-1] + dt*V;
float lambda = L*( sqrt(1 - (sin(theta)*(1-v*dt/L))**2) - cos(theta) )
// the array in which I stored the successive B points
Barray[i] = Barray[i-1] + lambda*u;
// the AB vector normalized
AiBiUnit = (Barray[i] - Aarray[i])/L;
// Refreshing the velocity of A
V = v*AiBiUnit;
// Refreshing u.
// u is indeed a unit vector separated from AiBiUnit by an angle theta,
// so you get it by rotating the newly computed AiBiUnit by an angle
// of +theta:
u = AiBiUnit.rotate(theta);
If you repeat a lot and/or increase the steering angle, the trajectory is a circle, which is coherent, I believe.