Okay, so this could be incredibly complicated or incredibly simple, depending on how one looks at it.
At the most basic level, you can treat this as a solution to the heat equation:
$$\frac{\partial T}{\partial t} = \alpha \nabla^2 T$$
where the Earth and the tropopause are the boundary conditions. In this model, we don't care about the transient, so the temporal derivative goes away. Since we don't care about the transient, we take the average temperature of the Earth over time (so the diurnal cycle-averaged temperature) and impose that value as the BC at $z=0$ and we know the beginning of the tropopause the temperature is again constant, so we set the BC at $z=11km$ to the temperature there.
This then gives you a linear solution for temperature, exactly what your lapse rate equation gives.
Now, your deviations in observations could be due to several things. First, the equation is for the time-averaged lapse rate; how many samples did you take? Second, the value of $\gamma$ depends heavily on humidity; what were the conditions at the times of the observations?
If you wanted to get fully into how to do all of this, it's not easy as others have alluded to. For regions under, roughly, 1200 m, we are within the atmospheric or planetary boundary layer. The composition of this layer is shown in the figure [from Wikipedia].
In this layer, there's a lot that goes on. But you can actually solve it in 2D (in spherical coordinates even) using the turbulent boundary layer equations, being sure of course to include the buoyancy terms and Coriolis terms if you are looking at a large enough scale where that matters. And really you need the low-Mach number form of them (rather than incompressible) to include the temperature equation and density variations.
Outside of 1200 meters or so, you can use the inviscid, low-Mach number Euler equations, again in 2D and again including buoyancy and Coriolis. You don't need to worry about the viscosity here to get decent answers, but know that from 1200 meters upwards is probably turbulent and that makes accurate solutions much, much harder. But for simplicity sake, you can keep it 2D and inviscid.
Then you need to blend the two solvers together at the edge of the boundary layer. The reason we do it this way is because the Euler equations+Boundary layer equations are much, much easier (less expensive) to solve even with the coupling of the two solutions through boundary conditions than it would be to get a fully turbulent Navier-Stokes solution everywhere.
The aerospace literature is very rich with how to do this and it would apply equally well to atmospheric work because the equations are the same (with the additional body force and source terms for rotating planets and buoyancy). A very quick internet search for "boundary layer coupled euler solver" turned up several recent AIAA papers from the same group, one of which is linked below for completeness.
A Rapid, Robust, and Accurate Coupled Boundary-Layer Method for Cart3D, AIAA-2012-0302
Edit:
I just re-read the question a bit more carefully and see that you are only taking measurements over the first few meters. I won't get full into boundary layer theory, but within the distance very close to the wall, nothing is linear at all. Here is an example of the profiles in a flat-plate boundary layer for velocity and temperature [from Wikipedia:
You can see that even on a perfectly flat plate with no roughness or variation in temperature along the plate, there is a non-linear region before the temperature becomes constant.
Again, using the boundary layer equations with the applicable terms in it will give you the solution you are looking for and it will certainly be non-linear near the surface, even for surfaces that are "simple" like an infinite, smooth flat plate. The surface roughness and non-uniform temperatures will only serve to make it less linear in that region.
