clc; clear all;
Lx=1.5; Ly=1; h=0.5; % space dimension & mesh size
T_left = 0; T_bottom = 0; % dirichet BC (left, bottom)
nx = Lx/h-1; ny = Ly/h; % number of node
T_right=zeros(Ly/h+1,1);
for i=1:length(T_right)
    T_right(i)=3*((i-1)*h)^3; % neumann BC (right)
end
Tn_upper=zeros(Lx/h+1,1);
for i=1:length(Tn_upper)
    Tn_upper(i)=6*((i-1)*h); % neumann BC (upper)
end
T_numbering = [];
iter = 0;
i = [1:nx]';
j = ones(nx,1);
for k = 1:ny
    T_numbering = [T_numbering;i,j*k];
end
f = zeros(nx*ny,1);
for j = 1:ny
    for i = 1:nx
        T_numbering(:,3) = zeros(nx*ny,1);
        a = find(T_numbering(:,1) == i+1 & T_numbering(:,2) == j);
        b = find(T_numbering(:,1) == i-1 & T_numbering(:,2) == j);
        c = find(T_numbering(:,1) == i & T_numbering(:,2) == j+1);
        d = find(T_numbering(:,1) == i & T_numbering(:,2) == j-1);
        e = find(T_numbering(:,1) == i & T_numbering(:,2) == j);
        T_numbering(a,3) = 1;
        T_numbering(b,3) = 1;
        T_numbering(c,3) = 1;
        T_numbering(d,3) = 1;
        T_numbering(e,3) = -4;
        iter = iter + 1;     
        T(iter,:) = T_numbering(:,3)';
        f(iter) = 12*(i*h)*(j*h)*h^2; % f value
        a = [a 1];
        b = [b 1];
        c = [c 1];
        d = [d 1];
        if length(a) == 1
            f(iter) = f(iter) - T_right(j+1);
        end
        if length(b) == 1
            f(iter) = f(iter) - T_left;
        end
        if length(c) == 1 % using temp. node 
            f(iter) = f(iter) - 2*h*Tn_upper(i+1); % calculation f
            T(iter,iter-nx)=T(iter,iter-nx)+1; % apply K matrix 
        end
        if length(d) == 1
            f(iter) = f(iter) - T_bottom;
        end
    end
end
Temp_temp1 = T\f;
iter = 0;
for i = ny:-1:1
    for j = 1:nx
        iter = iter + 1;
        Temp_temp2(i,j) = Temp_temp1(iter,1);
    end
end
Temp = zeros(ny+1,nx+2);
Temp(:,1) = T_left;
for i=1:ny+1
    Temp(i,nx+2) = T_right(ny+2-i);
end
Temp(ny+1,:) = T_bottom;
Temp(1:ny,2:nx+1) = Temp_temp2;
x=[0:h:h*(nx+1)]; y=[h*ny:-h:0];
surf(x,y,Temp)
xlabel('x'); ylabel('y');zlabel('Temp'); colorbar
% compare exact sol.(handbook) with numerical sol.
u11=Temp_temp1(find(T_numbering(:,1) == 0.5/h & T_numbering(:,2) == 0.5/h))
u21=Temp_temp1(find(T_numbering(:,1) == 1/h & T_numbering(:,2) == 0.5/h))
u12=Temp_temp1(find(T_numbering(:,1) == 0.5/h & T_numbering(:,2) == 1/h))
u22=Temp_temp1(find(T_numbering(:,1) == 1/h & T_numbering(:,2) == 1/h))