diff --git a/Project.toml b/Project.toml
index d437c3850..fe6660cb1 100644
--- a/Project.toml
+++ b/Project.toml
@@ -27,6 +27,7 @@ Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
diff --git a/demo/plotting/plot_inputs.jl b/demo/plotting/plot_inputs.jl
new file mode 100644
index 000000000..1168c0de9
--- /dev/null
+++ b/demo/plotting/plot_inputs.jl
@@ -0,0 +1,24 @@
+using UncertaintyQuantification, Plots
+
+X1 = RandomVariable(ProbabilityBox{Normal}(Dict(:μ => Interval(-1, 2), :σ => 1)), :X1)
+X2 = RandomVariable(ProbabilityBox{Normal}(Dict(:μ => Interval(-1, 2), :σ => 2)), :X2)
+X3 = IntervalVariable(-1, 2, :s)
+h = RandomVariable(Normal(0.24, 0.01), :h) # height
+
+plot(X1)    # p-box
+plot(X1, color = "red") # red p-box
+plot(X3)    # Interval
+plot(h)     # distribution 
+
+inputs = [X1, X2, X3, h]
+
+plot(inputs)    # Everything together
+
+samples = sample(inputs, 200)
+
+plot(X1)
+plot!(samples.X1)    # Samples ecdf samples of X1
+
+plot(samples.X1[1], samples.X2[1])  # Plots 2D box
+
+plot(samples.X1, samples.X2)        # Plots bivariate random sets of X1 and X2
diff --git a/src/UncertaintyQuantification.jl b/src/UncertaintyQuantification.jl
index 98d290e8c..073796fb0 100644
--- a/src/UncertaintyQuantification.jl
+++ b/src/UncertaintyQuantification.jl
@@ -21,6 +21,7 @@ using Random
 using Reexport
 using Roots
 using StatsBase
+using RecipesBase
 
 @reexport using Distributions
 
@@ -231,4 +232,11 @@ include("reliability/probabilityoffailure.jl")
 include("reliability/probabilityoffailure_imprecise.jl")
 include("sensitivity/sobolindices.jl")
 
+include("plotting/plot_recipes.jl")
+
+include("util/fourier-transform.jl")
+include("util/wrap.jl")
+include("util/imprecise.jl")
+include("util/kde.jl")
+
 end
diff --git a/src/inputs/imprecise/interval.jl b/src/inputs/imprecise/interval.jl
index 10ec31c34..88a146fdd 100644
--- a/src/inputs/imprecise/interval.jl
+++ b/src/inputs/imprecise/interval.jl
@@ -44,6 +44,11 @@ function bounds(i::Interval)
     return i.lb, i.ub
 end
 
+hi(X::Interval) = X.ub
+lo(X::Interval) = X.lb
+hi(X::Real) = X
+lo(X::Real) = X
+
 """
     IntervalVariable(lb::Real, ub::Real, name::Symbol)
 
diff --git a/src/plotting/plot_recipes.jl b/src/plotting/plot_recipes.jl
new file mode 100644
index 000000000..eeb3ac1d4
--- /dev/null
+++ b/src/plotting/plot_recipes.jl
@@ -0,0 +1,292 @@
+DEFAULT_ALPHA = 0.2
+DEFAULT_LABEL = ""
+DEFAULT_GRID = false
+DEFAULT_LEGEND = true
+DEFAULT_CDF = false
+DEFAULT_FILL_DISTRIBUTION=true
+DEFAULT_FILL_IMPRECISE=true
+DEFAULT_DISTRIBUTION = :pdf
+DEFAULT_FILL = :gray
+DEFAULT_COLOUR_PDF = :blue
+DEFAULT_COLOUR_UPPER = :red
+DEFAULT_COLOUR_LOWER = :black
+
+DEFAULT_INTERVAL_WIDTH=1.5
+DEFAULT_INTERVAL_EDGE_ALPHA=1
+
+DEFAULT_DISTRIBUTION_WIDTH=2
+
+DEFAULT_PLOT_RANGE_EXTEND_DENSITY = 0.2
+DEFAULT_PLOT_RANGE_EXTEND = 0.2
+DEFAULT_PLOT_RANGE_INTERVAL = 0.4
+DEFAULT_PLOT_GRID_NUMBER = 500
+DEFAULT_FONT_SIZE = 18
+DEFAULT_TICK_SIZE = 12
+
+###
+#   Plots for UQInputs
+###
+@recipe function _plot(
+    x::RandomVariable{T}; cdf_on=DEFAULT_CDF, shade=DEFAULT_FILL_DISTRIBUTION
+) where {T<:UnivariateDistribution}
+    grid --> DEFAULT_GRID
+    legend --> DEFAULT_LEGEND
+    label --> String(x.name)
+    (cdf_on ? ylabel --> "cdf" : ylabel --> "pdf")
+    seriescolor --> :auto
+
+    lo = quantile(x, 0.001)
+    hi = quantile(x, 0.999)
+    w = hi - lo
+    lo -= abs(w * DEFAULT_PLOT_RANGE_EXTEND_DENSITY)
+    hi += abs(w * DEFAULT_PLOT_RANGE_EXTEND_DENSITY)
+
+    xs = range(lo, hi, DEFAULT_PLOT_GRID_NUMBER)
+    ys = cdf_on ? cdf.(Ref(x), xs) : pdf.(Ref(x), xs)
+
+    # Primary line: either cdf or pdf
+    @series begin
+        seriestype := :path
+        fill := nothing
+        alpha --> 1
+        linewidth --> DEFAULT_DISTRIBUTION_WIDTH
+        label --> String(x.name)
+        xs, ys
+    end
+
+    # Optional fill for PDF only: reuse color, don't advance palette
+    if !cdf_on && shade
+        @series begin
+            primary := false          # <-- reuse the color from the primary line
+            seriestype := :path
+            fillrange := 0              # fill down to baseline
+            fillcolor := :match
+            fillalpha --> DEFAULT_ALPHA
+            linewidth --> 0             # fill only; no extra line
+            label := ""
+            xs, ys
+        end
+    end
+end
+
+@recipe function _plot(x::IntervalVariable; shade=DEFAULT_FILL_IMPRECISE)
+    # --- plot-level defaults (soft) ---
+    grid --> DEFAULT_GRID
+    legend --> DEFAULT_LEGEND
+    ylabel --> "cdf"
+    label --> String(x.name)             # single legend entry
+
+    seriescolor --> :auto
+
+    lo_grid = x.lb
+    hi_grid = x.ub
+
+    width = hi_grid - lo_grid
+
+    plot_lo = lo_grid - abs(width * DEFAULT_PLOT_RANGE_INTERVAL)    # For adding a slight width to the left and right side of the interval plot
+    plot_hi = hi_grid + abs(width * DEFAULT_PLOT_RANGE_INTERVAL)
+
+    x_grid = range(lo_grid, hi_grid, DEFAULT_PLOT_GRID_NUMBER)
+
+    cdf_lo = x_grid .>= x.ub
+    cdf_hi = x_grid .> x.lb
+
+    # Plot upper cdf (primary, inherits colour and label)
+    @series begin
+        seriestype := :path
+        alpha --> 1
+        linewidth --> DEFAULT_DISTRIBUTION_WIDTH
+        x_grid, cdf_hi
+    end
+
+    # Plot lower cdf
+    @series begin
+        primary := false
+        seriestype := :path
+        alpha --> 1
+        label := ""
+        linewidth --> DEFAULT_DISTRIBUTION_WIDTH
+        x_grid, cdf_lo
+    end
+
+    # Add invisible lines so interval bounds don't touch plot boundaries
+    @series begin
+        primary := false
+        seriestype := :path
+        alpha --> 0
+        label := ""
+        linewidth --> 0
+        range(plot_lo, plot_hi, 100), zeros(100)
+    end
+
+    # Plot fill
+    if shade
+        @series begin
+            primary := false
+            seriestype := :path
+            fillcolor := :match              # match this series' line color
+            fillrange := cdf_hi
+            color := DEFAULT_FILL
+            fillalpha := DEFAULT_ALPHA
+            linewidth --> 0                  # draw fill only here
+            label := ""
+            x_grid, cdf_lo
+        end
+    end
+end
+
+@recipe function _plot(
+    x::RandomVariable{T}; shade=DEFAULT_FILL_IMPRECISE
+) where {T<:ProbabilityBox}
+    # --- plot-level defaults (soft) ---
+    grid --> DEFAULT_GRID
+    legend --> DEFAULT_LEGEND
+    ylabel --> "cdf"
+    label --> String(x.name)             # single legend entry
+
+    seriescolor --> :auto
+
+    lo_grid = quantile(x, 0.001).lb
+    hi_grid = quantile(x, 0.999).ub
+    width = hi_grid - lo_grid
+    lo_grid = lo_grid - abs(width * DEFAULT_PLOT_RANGE_EXTEND)
+    hi_grid = hi_grid + abs(width * DEFAULT_PLOT_RANGE_EXTEND)
+
+    x_grid = range(lo_grid, hi_grid, DEFAULT_PLOT_GRID_NUMBER)
+    cdf_evals = cdf.(Ref(x), x_grid)
+
+    # Plot upper cdf (primary, inherits colour and label)
+    @series begin
+        seriestype := :path
+        alpha --> 1
+        linewidth --> DEFAULT_DISTRIBUTION_WIDTH
+        x_grid, hi.(cdf_evals)
+    end
+
+    # Plot lower cdf
+    @series begin
+        primary := false
+        seriestype := :path
+        alpha --> 1
+        linewidth --> DEFAULT_DISTRIBUTION_WIDTH
+        label := ""
+        x_grid, lo.(cdf_evals)
+    end
+
+    # Plot fill
+    if shade
+        @series begin
+            primary := false
+            seriestype := :path
+            fillcolor := :match              # match this series' line color
+            fillrange := hi.(cdf_evals)
+            fillalpha --> DEFAULT_ALPHA
+            linewidth --> 0                  # draw fill only here
+            label := ""
+            x_grid, lo.(cdf_evals)
+        end
+    end
+end
+
+@recipe function _plot(x::Vector{T}) where {T<:UQInput}
+    # Filter out Parameter objects
+    x_no_params = filter(xi -> !isa(xi, Parameter), x)
+
+    N = length(x_no_params)
+    cols = ceil(Int, sqrt(N))
+    rows = ceil(Int, N / cols)
+    layout --> (rows, cols)
+
+    # Choose a grid palette once (users can still override via plot(...; palette=...))
+    # palette --> :default  # or :default, :Dark2_8, etc.
+
+    for i in 1:N
+        @series begin
+            subplot := i
+            seriescolor --> i  # <-- panel i uses the i-th color from the current palette
+            x_no_params[i]
+        end
+    end
+end
+
+###
+#   This code is a modified version of the plot recipe from IntervalArithmetic.jl
+#       https://github.com/JuliaIntervals/IntervalArithmetic.jl       
+###
+
+# Plot a 2D IntervalBox:
+@recipe function _plot(x::Interval, y::Interval)
+    seriesalpha --> DEFAULT_ALPHA
+    seriestype := :shape
+
+    label := false
+
+    linecolor --> :black                        # Explicitly set edge color
+    linewidth --> DEFAULT_INTERVAL_WIDTH        # Make edges more visible
+    linealpha --> DEFAULT_INTERVAL_EDGE_ALPHA
+
+    x = [x.lb, x.ub, x.ub, x.lb]
+    y = [y.lb, y.lb, y.ub, y.ub]
+
+    x, y
+end
+
+# Plot a vector of 2D IntervalBoxes:
+@recipe function _plot(xx::Vector{T}, yy::Vector{T}) where {T<:Interval}
+    seriesalpha --> DEFAULT_ALPHA
+    seriestype := :shape
+
+    label := false
+
+    linecolor := :black                         # Explicitly set edge color
+    linewidth --> DEFAULT_INTERVAL_WIDTH        # Make edges more visible
+    linealpha --> DEFAULT_INTERVAL_EDGE_ALPHA
+
+    xs = Float64[]
+    ys = Float64[]
+
+    # build up coords:  # (alternative: use @series)
+    for i in 1:length(xx)
+        (x, y) = (xx[i], yy[i])
+
+        # use NaNs to separate
+        append!(xs, [x.lb, x.ub, x.ub, x.lb, NaN])
+        append!(ys, [y.lb, y.lb, y.ub, y.ub, NaN])
+    end
+
+    xs, ys
+end
+
+###
+# Plots for samples of data frames
+###
+
+@recipe function _plot(x::Vector{Interval})
+    if length(unique(x))==1
+        return x[1]
+    else
+        grid --> DEFAULT_GRID
+        legend --> DEFAULT_LEGEND
+
+        # xlabel --> x[1].name
+        ylabel --> "cdf"
+        N_samples = length(x)
+
+        lows = sort(lo.(x))
+        his = sort(hi.(x))
+
+        is = range(0, 1, length=N_samples)
+
+        @series begin
+            seriestype := :steppre
+            color --> DEFAULT_COLOUR_LOWER
+            lows, is
+        end
+
+        @series begin
+            seriestype := :steppost
+            color --> DEFAULT_COLOUR_UPPER
+            his, is
+        end
+    end
+end
diff --git a/test/Project.toml b/test/Project.toml
index f36f601fe..af12538e8 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -4,6 +4,7 @@ Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 HCubature = "19dc6840-f33b-545b-b366-655c7e3ffd49"
 HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
diff --git a/test/inputs/imprecise/interval.jl b/test/inputs/imprecise/interval.jl
index 54d202a10..6c0fe10a9 100644
--- a/test/inputs/imprecise/interval.jl
+++ b/test/inputs/imprecise/interval.jl
@@ -1,3 +1,5 @@
+using UncertaintyQuantification: hi, lo
+
 @testset "Interval" begin
     name = :l
     lb = 0.14
@@ -16,6 +18,12 @@
     @test !(0.17 ∈ interval)
 
     @test sprint(show, interval) == "[0.14, 0.16]"
+
+    @test hi(interval) == interval.ub
+    @test lo(interval) == interval.lb
+
+    @test hi(2.0) == 2.0
+    @test lo(2.0) == 2.0
 end
 
 @testset "IntervalVariable" begin
diff --git a/test/plotting/plotting.jl b/test/plotting/plotting.jl
new file mode 100644
index 000000000..200743349
--- /dev/null
+++ b/test/plotting/plotting.jl
@@ -0,0 +1,59 @@
+@testset "Plotting Recipes" begin
+
+    @testset "RandomVariable PDF Plot" begin
+        rv = RandomVariable(Normal(0, 1), :X)
+        plt = plot(rv; cdf_on=false)
+        @test typeof(plt) <: Plots.Plot
+        @test plt.series_list[1][:seriestype] == :path
+        @test plt.series_list[2][:seriestype] == :path
+    end
+
+    @testset "RandomVariable CDF Plot" begin
+        rv = RandomVariable(Normal(0, 1), :X)
+        plt = plot(rv; cdf_on=true)
+        @test typeof(plt) <: Plots.Plot
+    end
+
+    @testset "IntervalVariable Plot" begin
+        iv = IntervalVariable(1.0, 2.0, :Y)
+        plt = plot(iv)
+        @test typeof(plt) <: Plots.Plot
+        @test length(plt.series_list) == 3
+    end
+
+    @testset "ProbabilityBox Plot" begin
+        pb = RandomVariable(ProbabilityBox{Normal}(Dict(:μ => Interval(-1, 2), :σ => 2)), :X2)
+        plt = plot(pb)
+        @test typeof(plt) <: Plots.Plot
+        @test length(plt.series_list) == 3
+    end
+
+    @testset "Vector of UQInputs Plot" begin
+        inputs = [RandomVariable(Normal(0,1), :A), IntervalVariable(1.0, 2.0, :B)]
+        plt = plot(inputs)
+        @test typeof(plt) <: Plots.Plot
+        @test plt.layout isa Plots.GridLayout
+    end
+
+    @testset "2D IntervalBox Plot" begin
+        x = Interval(1.0, 2.0)
+        y = Interval(3.0, 4.0)
+        plt = plot(x, y)
+        @test typeof(plt) <: Plots.Plot
+        @test plt.series_list[1][:seriestype] == :shape
+    end
+
+    @testset "Vector of IntervalBoxes Plot" begin
+        xs = [Interval(1.0, 2.0), Interval(2.0, 3.0)]
+        ys = [Interval(3.0, 4.0), Interval(4.0, 5.0)]
+        plt = plot(xs, ys)
+        @test typeof(plt) <: Plots.Plot
+        @test plt.series_list[1][:seriestype] == :shape
+    end
+
+    @testset "Vector of Intervals Plot" begin
+        intervals = [Interval(1.0, 2.0), Interval(1.5, 2.5), Interval(2.0, 3.0)]
+        plt = plot(intervals)
+        @test typeof(plt) <: Plots.Plot
+    end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index eb3e174ef..499938f97 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -7,6 +7,7 @@ using QuasiMonteCarlo
 using Random
 using StatsBase: fit, Histogram, corkendall
 using Test
+using Plots
 using UncertaintyQuantification
 
 include("inputs/empiricaldistribution.jl")
@@ -49,6 +50,8 @@ include("simulations/subset.jl")
 
 include("util/fourier-transform.jl")
 
+include("plotting/plotting.jl")
+
 if Sys.islinux()
     HPC = false
     HPC_account = "HPC_account_1"