diff --git a/R/internal.R b/R/internal.R
index 21b2c8d..e4a593e 100644
--- a/R/internal.R
+++ b/R/internal.R
@@ -25,18 +25,12 @@ set_xout <- function(chartcode, yname) {
   if (design == "A") {
     return(round(seq(0.5, 15, 0.5) / 12, 4L))
   }
-  if (design == "B" & yname == "wfh") {
-    return(round(seq(50, 120, by = 2), 4L))
-  }
   if (design == "B" & yname %in% c("hgt", "dsc")) {
     return(round(c(0.5, 0.75, 1:48) / 12, 4L))
   }
-  if (design == "B" & yname == "hdc") {
+  if (design == "B" & yname %in% c("hdc", "wfh")) {
     return(round(seq(0.1, 4, by = 0.1), 4L))
   }
-  if (design == "C" & yname == "wfh") {
-    return(round(seq(60, 184, by = 4), 4L))
-  }
   if (design == "C") {
     return(round(seq(1, 21, by = 0.5), 4L))
   }
diff --git a/R/set_curves.R b/R/set_curves.R
index 00f40f5..f6bc245 100644
--- a/R/set_curves.R
+++ b/R/set_curves.R
@@ -37,9 +37,13 @@ set_curves <- function(g,
     mutate(
       id = -1,
       sex = child$sex,
-      ga = child$ga
+      ga = child$ga,
+      x2 = .data$x
     ) %>%
-    select(all_of(c("id", "age", "xname", "yname", "x", "y", "sex", "ga")))
+    select(all_of(c("id", "age", "xname", "yname", "x", "y", "sex", "ga", "x2")))
+  # For WFH, temporary sort on age to get correct time sequence, use x2 to store x
+  idx <- data$yname == "wfh"
+  data$x[idx] <- data$age[idx]
 
   # get data of matches
   time <- vector("list", length(ynames))
@@ -107,7 +111,17 @@ set_curves <- function(g,
   # append synthetic data
   data <- data %>%
     bind_rows(synt) %>%
-    arrange(.data$id, .data$yname, .data$x, .data$age) %>%
+    arrange(.data$id, .data$yname, .data$x)
+
+  # For wfh, interpolate hgt from age, and overwrite data$x with hgt
+  idx <- data$yname == "wfh"
+  if (any(idx)) {
+    data$x[idx] <- safe_approx(x = data$x[idx], y = data$x2[idx],
+                               xout = data$x[idx], ties = list("ordered", mean))$y
+  }
+
+  # add Z-scores
+  data <- data %>%
     select(-"age") %>%  # fool set_refcodes()
     mutate(refcode_z = nlreferences::set_refcodes(.)) %>%
     mutate(
@@ -133,8 +147,7 @@ set_curves <- function(g,
     group_by(.data$id, .data$yname, .data$pred) %>%
     mutate(z = safe_approx(
       x = .data$x, y = .data$z, xout = .data$x,
-      ties = list("ordered", mean)
-    )$y) %>%
+      ties = mean)$y) %>%
     ungroup()
 
   # set refcode as target's sex and ga
diff --git a/tests/testthat/test-process_chart.R b/tests/testthat/test-process_chart.R
index f941c8d..310491a 100644
--- a/tests/testthat/test-process_chart.R
+++ b/tests/testthat/test-process_chart.R
@@ -22,6 +22,28 @@ ind <- bdsreader::read_bds(fn)
 g <- process_chart(ind, chartcode = "NJAA",
                    dnr = "smocc", period = c(0.5, 1.1667), nmatch = 10, break_ties = TRUE)
 
+# incorrect order of observations for WFH when hgt is not monotone
+day <- c(0, 13, 42, 91, 152, 287, 336, 434, 541, 632, 744, 905)
+hgt <- c(NA, NA, 56, 61.5, 67, 72.5, 74, 78, 83, 84, 89, 88)
+wgt <- c(2.879, 3.14, 4.4, 6.055, 7.15, 7.915, 8.25, 9.45, 10.8, 10.45, 10.8, 11.9)
+df <- data.frame(age = round(day / 365.25, 4), hgt, wgt)
+xyz <- ind$xyz[ind$xyz$yname == "wfh", ]
+xyz$age <- df$age[-(1:2)]
+xyz$x <- df$hgt[-(1:2)]
+xyz$y <- df$wgt[-(1:2)]
+xyz$z <- centile::y2z(y = xyz$y, x = xyz$x, refcode = "nl_1997_wfh_female_nla",
+                      pkg = "nlreferences", rule = 2L)
+data <- ind
+data$xyz <- xyz
+
+# Note WFH curve age sequence: correct for curve_interpolation is FALSE
+# But has approximation errors when curve_interpolation is TRUE
+test_that("Weight for height curve has correct time sequence", {
+  expect_silent(process_chart(data, chartcode = "NJBR",
+                              dnr = "smocc", period = c(0.5, 1.1667), nmatch = 10,
+                              break_ties = TRUE, curve_interpolation = FALSE))
+})
+
 test_that("terneuzen donordata yields matches", {
   expect_silent(process_chart(ind, chartcode = "NJCH",
                               dnr = "terneuzen", period = c(2, 18),