Procedural heightmap shader - Normals and tangents

Hi, I've been working on a procedural terrain shader and though the heightmap generates correctly, I'm having trouble with shading.

Let's assume I have the following function : float height = f(position). This function returns a height value (y) based on the given position (x,z). My shader looks like this :

 struct Input 
 {
     float2 uv_MainTex;
     float3 worldPos;
 };

 void vert (inout appdata_full v, out Input o) 
 {
     o.worldPos = mul(_Object2World, v.vertex).xyz;
     float2 position = o.worldPos.xz;
     float height = f(position);
     v.vertex.y = height;
     o.worldPos = mul(_Object2World, v.vertex).xyz;
 }

 void surf (Input IN, inout SurfaceOutput o) 
 {
     o.Albedo = tex2D(_MainTex, IN.uv_MainTex).rgb;
 }

I'm using the worldPos so that the terrain is continuous over several grids (each is 192x192 triangle pairs). This shader outputs the terrain correctly, but it is not shaded as you can imagine.

This is where I have a problem : how can I calculate normals/tangents for the generated surface? Note that an analytical approach using derivatives is NOT an option. I tried a numerical approach by sampling the heights in neighbouring positions but it's giving me weird results and doesn't compute tangents. This is the code I used (in the surf function) :

 // Compute Normals
 float d = 0.1;
 float2 grad;
 float2 position = IN.worldPos;
 grad.x = f(float2(position.x+d, position.y)) -
          f(float2(position.x-d, position.y));
 grad.y = f(float2(position.x, position.y+d)) -
          f(float2(position.x, position.y-d));
 o.Normal = normalize(cross(float3(1.0, 0.0, grad.x), float3(0.0, 1.0, grad.y)));

Where d is a small offset to the actual position. I know there is something wrong with this code but I'm out of ideas.

Any help would be greatly appreciated!