using Plots, LinearAlgebra, SparseArrays

uexact(x,y)=y/((1+x)^2+y^2)

M=100
h=1/(M+1)

exact=zeros(M+2,M+2)
for k=0:M+1
    for l=0:M+1
        exact[k+1,l+1]=uexact(k*h,l*h)
    end
end

xs=(0:M+1)*h

#heatmap(xs,xs,exact')
#contour(xs,xs,exact')
#surface(xs,xs,exact')

function phi(x,y)
    if y==0
        return 0.0
    elseif x==0
        return y/(1+y^2)
    elseif y==1
        return 1/((1+x)^2+1)
    elseif x==1
        return y/(4+y^2)
    end
    println("the point ($x,$y) was not on the boundary!")
end

function f(x,y)
    return 0
end

# Mapping between the grid points and how they are
# stored:  u(x_k,y_l)=u[l+(k-1)*M]
function klton(k,l)
    return l+(k-1)*M
end

M2=M^2
A=sparse(1:M2,1:M2,-4*ones(M2))
c=zeros(M2)

for k=1:M
    for l=1:M
        c[klton(k,l)]=h^2*f(k*h,l*h)
    end
end

function interior(k,l)
    if k>=1&&k<=M&&l>=1&&l<=M
        return true
    end
    return false
end

for k=1:M
    for l=1:M
        function update(k1,l1)
            if interior(k1,l1)
                A[klton(k,l),klton(k1,l1)]=1
            else
                c[klton(k,l)]-=phi(k1*h,l1*h)
            end
        end
        update(k,l+1)            
        update(k,l-1)
        update(k+1,l)
        update(k-1,l)
    end
end

u=A\c

approx=zeros(M+2,M+2)
for k=0:M+1
    for l=0:M+1
        if interior(k,l)
            approx[k+1,l+1]=u[klton(k,l)]
        else
            approx[k+1,l+1]=phi(k*h,l*h)
        end
    end
end

#heatmap(xs,xs,approx')
#contour(xs,xs,approx')
#surface(xs,xs,approx')

heatmap(xs,xs,1000*(approx-exact)')