Solving heat equation with Dirichlet boundary

Hey, I'm solving the heat equation on a grid for time with inhomogeneous Dirichlet boundary conditions . I'm using the implicit scheme for FDM, so I'm solving the Laplacian with the five-point-stencil, i.e. where are indices of the mesh. With the implicit scheme for the heat equation we get to solve where A is the matrix representing the discretized Laplacian, and F is zero if is in the middle of the grid and on the boundary because of the Dirichlet boundary which needs to be included in the scheme. My problem is, that the boundary of my plot is still equal to 0 and not and I don't know why since I've changed the right hand side F according to the boundary. Even if I initialize U with instead of I get the boundary 0. Does somebody know why this happens?

Here is my code:

function diffusion
%heat equation is solved with lin. system of equations
T = 0.01;
J = 200; K = 200; L = (J-1)^2;
Delta_x = 1/J; Delta_t = 1/K;   %mesh-size, timestep-size
lambda = Delta_t/Delta_x^2;     %
end_time = floor(T/Delta_t);
R0 = 0.3;                       %radius of the circle
U = zeros(K+1,L);               
e = ones(J-1,1); E = ones(L,1);
X = spdiags([-e,4*e,-e],[-1,0,1],J-1,J-1);  %matrix to the left and right point in 5-point-stencil
A = sparse(L,L);                            %matrix for whole 5-point-stencil, it is a matrix which contains matrices
F = zeros(1,L);                             %right hand side
for j=1:J-2
    F((J-1)*j) = -1;                        %points on boundary of the grid, but not corner
    F((J-1)*j+1) = -1;
end
F(1) = -2; F(J-1) = -2; F((J-2)*(J-1)+1) = -2; F(L) = -2;   %points on each corner of the grid
for j = 1:(J-1):L
    A(j:j+J-2,j:j+J-2) = X;                 %puts X on the diagonal of A
end
A = A+spdiags([-E,-E],[-J+1,J-1],L,L);      %left and right we get -E because of upper and downer points of 5-point-stencil
B = speye(L)+lambda*A;
C = speye(L);
%initial condition
for j = 1:J-1
    for m =1:J-1
        U(1,j+(m-1)*(J-1)) = u_0([j,m]*Delta_x,R0);
    end
end

for k = 0:end_time
    show(U(k+1,:),J,Delta_x); drawnow;
    U(k+2,:) = B\(C*U(k+1,:)'+lambda*F');       %solve the implicit scheme
end

function val = u_0(x,R)
%initial condition is a circle
if ((x(1)-1/2)^2+(x(2)-1/2)^2)<=R^2
    u = 1;
else
    u = -1;
end
val = u;

function show(U,J,Delta_x)
U_mat = zeros(J+1,J+1); U_mat(2:J,2:J) = reshape(U,J-1,J-1)';
mesh(Delta_x*(0:J),Delta_x*(0:J),U_mat); %axis([0,1,0,1,0,1,0,1]); 

Here is the plot with the boundary still at 0 instead of .

[darova]

It's too complicated. I can only attach some example, maybe you will find it helpfull