The equation given, (d gamma / d u)|_{u=u0} = tau, doesn't define a differential equation for gamma, it defines a boundary condition: gamma dot = tau at u=u0.
If you mean you've got discrete lists of [tau_i gamma_i], Matlab's interp1 will create the appearance of a smooth function that can then be used with e.g. fzero to search for the zero.
Matlab's ODE solvers will give output for specific times if you specify them,
[tout, yout] = ode*(odefun, [t0 t1 t2 ... tfinal], y0)
where the vector of t_i are monotonically increasing or decreasing (see doc ode45 for more)
