diff --git a/docs/pages.jl b/docs/pages.jl
index 24d4d3a7..455a9702 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -15,6 +15,7 @@ pages = ["index.md"
"Variable Fidelity" => "variablefidelity.md",
"Gradient Enhanced Kriging" => "gek.md",
"GEKPLS" => "gekpls.md",
+ "Improved GEKPLS" => "Improvedgekpls.md",
"MOE" => "moe.md",
"Parallel Optimization" => "parallel.md",
]
diff --git a/docs/src/Improvedgekpls.md b/docs/src/Improvedgekpls.md
new file mode 100644
index 00000000..52ac943e
--- /dev/null
+++ b/docs/src/Improvedgekpls.md
@@ -0,0 +1,55 @@
+# GEKPLS Function
+
+Gradient Enhanced Kriging with Partial Least Squares Method (GEKPLS) is a surrogate modelling technique that brings down computation time and returns improved accuracy for high-dimensional problems. The Julia implementation of GEKPLS is adapted from the Python version by [SMT](https://github.com/SMTorg) which is based on this [paper](https://arxiv.org/pdf/1708.02663.pdf).
+
+# Modifications for Improved GEKPLS Function:
+
+To enhance the GEKPLS function, sampling method was changed from ```SobolSample()``` to ```HaltonSample()```.
+
+
+```@example gekpls_water_flow
+
+using Surrogates
+using Zygote
+
+function water_flow(x)
+ r_w = x[1]
+ r = x[2]
+ T_u = x[3]
+ H_u = x[4]
+ T_l = x[5]
+ H_l = x[6]
+ L = x[7]
+ K_w = x[8]
+ log_val = log(r/r_w)
+ return (2*pi*T_u*(H_u - H_l))/ ( log_val*(1 + (2*L*T_u/(log_val*r_w^2*K_w)) + T_u/T_l))
+end
+
+n = 1000
+lb = [0.05,100,63070,990,63.1,700,1120,9855]
+ub = [0.15,50000,115600,1110,116,820,1680,12045]
+x = sample(n,lb,ub,HaltonSample())
+grads = gradient.(water_flow, x)
+y = water_flow.(x)
+n_test = 100
+x_test = sample(n_test,lb,ub,GoldenSample())
+y_true = water_flow.(x_test)
+n_comp = 2
+delta_x = 0.0001
+extra_points = 2
+initial_theta = [0.01 for i in 1:n_comp]
+g = GEKPLS(x, y, grads, n_comp, delta_x, lb, ub, extra_points, initial_theta)
+y_pred = g.(x_test)
+rmse = sqrt(sum(((y_pred - y_true).^2)/n_test)) #root mean squared error
+println(rmse) #0.0347
+```
+
+
+
+
+
+
+| **Sampling Method** | **RMSE** | **Differences** |
+|----------------------|--------------------------|---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
+| **Sobol Sampling** | 0.021472963465423097 | Utilizes digital nets to generate quasi-random numbers, offering low discrepancy points for improved coverage. - Requires careful handling, especially in higher dimensions. |
+| **Halton Sampling** | 0.02144270998045834 | Uses a deterministic sequence based on prime numbers to generate points, allowing for quasi-random, low-discrepancy sampling. - Simpler to implement but may exhibit correlations in some dimensions affecting coverage.
diff --git a/docs/src/tensor_prod.md b/docs/src/tensor_prod.md
index 832603cc..504c7a78 100644
--- a/docs/src/tensor_prod.md
+++ b/docs/src/tensor_prod.md
@@ -1,6 +1,9 @@
# Tensor product function
-The tensor product function is defined as:
-``f(x) = \prod_{i=1}^d \cos(a\pi x_i)``
+A tensor product function combines multiple functions or vectors using the tensor product operation. The tensor product is a mathematical operation that takes two vectors and produces another vector space, capturing their joint behavior across multiple dimensions.
+
+For instance, consider a tensor product function defined as follows:
+
+```\[ f(x) = ∏ᵢ=₁ᵈ cos(aπxᵢ) \]```
Let's import Surrogates and Plots:
```@example tensor
@@ -38,3 +41,15 @@ plot!(xs,f.(xs), label="True function", legend=:top)
plot!(xs, loba_1.(xs), label="Lobachevsky", legend=:top)
plot!(xs, krig.(xs), label="Kriging", legend=:top)
```
+
+## Kriging Plot
+
+
+
+## Lobachevsky Plot
+
+
+
+## Combined Plot
+
+