diff --git a/Project.toml b/Project.toml
index 516f264..57f1d8e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -3,6 +3,11 @@ uuid = "ce1ef771-b552-5ca6-80e4-7358b8c578e2"
 authors = ["Florian Oswald <florian.oswald@gmail.com>"]
 version = "0.1.0"
 
+[deps]
+ApproXD = "f0e97480-066f-5960-9ba5-e027352bc8af"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
diff --git a/q1.png b/q1.png
new file mode 100644
index 0000000..f3978bb
Binary files /dev/null and b/q1.png differ
diff --git a/src/HWfuncapp.jl b/src/HWfuncapp.jl
index 732fe3f..83d6576 100644
--- a/src/HWfuncapp.jl
+++ b/src/HWfuncapp.jl
@@ -4,21 +4,60 @@ using FastGaussQuadrature # to get chebyshevnodes
 
 # you dont' have to use PyPlot. I did much of it in PyPlot, hence
 # you will see me often qualify code with `PyPlot.plot` or so.
-using PyPlot  
+using PyPlot
 import ApproXD: getBasis, BSpline
 using Distributions
 using ApproxFun
 using Plots
 
+
 ChebyT(x,deg) = cos(acos(x)*deg)
 unitmap(x,lb,ub) = 2 .* (x .- lb) ./ (ub .- lb) .- 1	#[a,b] -> [-1,1]
 
 function q1(n)
 	f(x) = x .+ 2x.^2 - exp.(-x)
+# Generate Chebychev nodes
+	deg = n:-1:1
+	z = 3*cos.((2 .* deg .- 1).*pi ./ 2n)
+
+	# Compute function values at interpolation nodes
+	y = zeros(n)
+	for i in 1:n
+       y[i] = f(z[i])
+    end
+
+	# Compute basis matrix at the nodes
+	P = zeros(n, n)
+	for i in 0:n-1
+	       for j in 1:n
+	       P[j,i+1] = ChebyT.(unitmap(z, -3, 3)[j], i)
+	       end
+	end
+	# Obtain coefficients (solve the system)
+	c = inv(P) * y
+
+	# Recompute matrix for new set of points
+	# Evaluate Φ(x'): i.e. evaluate all 15 Chebyshev polynomials at the new x, and multiply them by c
+	n_new = 100
+	b = unitmap(range(-3, stop = 3, length = n_new), -3, 3)
+	P_b = zeros(n_new, n)
+	for i in 0:n-1
+	       for j in 1:n_new
+	       P_b[j,i+1] = ChebyT(b[j], i)
+	       end
+	end
+
+	# Multiply this new matrix by c to obtain approximation of f
+		y_new = f.(range(-3, stop = 3, length = n_new))
+		y_approx = (P_b) * c
+
+		approx_error = y_new .- y_approx
+		plot(range(-3, stop = 3, length = n_new), [y_new approx_error], label = ["Actuals" "Errors"], layout = 2)
+		scatter!(range(-3, stop = 3, length = n_new), y_approx, label = "Approximation", marker = (1, :circle))
 
 	# without using PyPlot, just erase the `PyPlot.` part
-	PyPlot.savefig(joinpath(dirname(@__FILE__),"..","q1.png"))
-	return Dict(:error=>maximum(abs,err))
+	savefig(joinpath(dirname(@__FILE__),"..","q1.png"))
+	return Dict(:error=>maximum(abs,approx_error))
 end
 
 function q2(b::Number)
@@ -43,7 +82,7 @@ end
 
 # I found having those useful for q5
 mutable struct ChebyType
-	f::Function # fuction to approximate 
+	f::Function # fuction to approximate
 	nodes::Union{Vector,LinRange} # evaluation points
 	basis::Matrix # basis evaluated at nodes
 	coefs::Vector # estimated coefficients
@@ -79,7 +118,7 @@ function q5(deg=(5,9,15),lb=-1.0,ub=1.0)
 
 	runge(x) = 1.0 ./ (1 .+ 25 .* x.^2)
 
-	
+
 	PyPlot.savefig(joinpath(dirname(@__FILE__),"..","q5.png"))
 
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 6453446..d635c7f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -3,7 +3,8 @@ using Test
 
 @testset "HWfuncapp.jl" begin
 	# test that q1 with 15 nodes has an error smaller than 1e-9
-    
+	a = q1(15)
+	@test a[:error] < 1e-9
     # test that the integral of function h [-10,0] ≈ 1.46039789878568
 
     @test false  # by default fail