# L-34 MCS 471 Wed 9 Nov 2022 : linearbvp2.jl

# Solving a linear boundary value problem with finite differences:
# y" = y, for x in [1, 3], y'(1) = 1.17520 and y'(3) = 10.0179.

using Printf
using LinearAlgebra

Base.show(io::IO, f::Float64) = @printf(io, "%.3e", f)

"""
    A, b = setup(dim::Int)

returns the linear system to solve a BVP 
with finite differences over the interval [a,b]
with dim points, with a step size h of (b-a)/(dim-1).
"""
function setup(dim::Int, a::Float64, b::Float64)
    h = (b - a)/(dim-1)
    d = -(2 + h^2)*ones(dim)
    e = ones(dim-1); e[1] = 2
    f = ones(dim-1); f[dim-1] = 2
    A = diagm(d) + diagm(+1 => e) + diagm(-1 => f)
    b = zeros(dim)
    yprime1 = 1.17520
    yprime3 = 10.0179
    b[1] = 2*h*yprime1
    b[dim] = -2*h*yprime3
    return A, b
end

"""
    Runs a particular case for dim points.
"""
function run(dim::Int, verbose::Bool=true)
    A, b = setup(dim, 1.0, 3.0)
    if verbose
        println("The matrix A :")
        show(stdout, "text/plain", A); println("")
        println("The vector b :")
        show(stdout, "text/plain", b); println("")
    end
    sol = A\b
    if verbose
        println("The solution :")
        show(stdout, "text/plain", sol); println("")
    end
    h = 2.0/(dim-1)
    sdim = @sprintf("%3d", dim)
    print("dim = ", sdim, " left slope : ")
    show(stdout, "text/plain", (sol[2] - sol[1])/h);
    print(", right slope : ")
    show(stdout, "text/plain", (sol[dim] - sol[dim-1])/h); println("")
end

function main()
   run(5)
   run(20,false)
   run(100,false)
   run(800,false)
end

main()