function w=poisson_inputpower_half(xl,xr,yb,yt,M,N,P,K)
H=0.005;
% K=3.85;%$for step4
delta=0.1;
u_b=20;
L=2;
% f=@(x,y) 0; % define input function data 
% g1=@(x) 0; %bottom % define boundary values 
% g2=@(x) 0; %top % Example 8.8 is shown 
% g3=@(y) 0; %left
% g4=@(y) 0; %right
m=M+1;n=N+1; mn=m*n; 
h=(xr-xl)/M;h2=h^2;k=(yt-yb)/N;k2=k^2;
x=xl+(0:M)*h; % set mesh values 
y=yb+(0:N)*k; 
A=zeros(mn,mn);b=zeros(mn,1); 
for i=2:m-1 % interior points 
  for j=2:n-1 
    A(i+(j-1)*m,i-1+(j-1)*m)=1/h2;
    A(i+(j-1)*m,i+1+(j-1)*m)=1/h2; 
    A(i+(j-1)*m,i+(j-1)*m)=-2/h2-2/k2-2*H/(K*delta); 
    A(i+(j-1)*m,i+(j-2)*m)=1/k2;
    A(i+(j-1)*m,i+j*m)=1/k2; 
    b(i+(j-1)*m)=0; 
  end 
end 
for i=1:m % bottom and top boundary points 
  j=1;
  A(i+(j-1)*m,i+(j-1)*m)=-3/(2*k)+H/K;
  A(i+(j-1)*m,i+(j)*m)=4/(2*k);
  A(i+(j-1)*m,i+(j+1)*m)=-1/(2*k);
  b(i+(j-1)*m)=0; 
%   j=2;
%   A(i+(j-1)*m,i+(j)*m)=4/(2*k);
%   b(i+(j-1)*m)=0; 
%   j=3;
%   A(i+(j-1)*m,i+(j+1)*m)=-1/(2*k);
%   b(i+(j-1)*m)=0; 
%   j=n-2;
%   A(i+(j-1)*m,i+(j-3)*m)=-1/(2*k);
%   b(i+(j-1)*m)=0; 
%   j=n-1;
%   A(i+(j-1)*m,i+(j-2)*m)=4/(2*k);
%   b(i+(j-1)*m)=0; 
  j=n;
  A(i+(j-1)*m,i+(j-1)*m)=-3/(2*k)+H/K;
  A(i+(j-1)*m,i+(j-2)*m)=4/(2*k);
  A(i+(j-1)*m,i+(j-3)*m)=-1/(2*k);
  b(i+(j-1)*m)=0; 
end 
%power input
for j=2:(n-1)/2 % left and right boundary points 
%   i=1;
%   A(i+(j-1)*m,i+(j-1)*m)=-3/(2*h)+H/K;
%   b(i+(j-1)*m)=g3(y(j)); 
%   i=2;
%   A(i+(j-1)*m,i+(j-1)*m)=4/(2*h);
%   b(i+(j-1)*m)=g3(y(j));
%   i=3;
%   A(i+(j-1)*m,i+(j-1)*m)=-1/(2*h);
%   b(i+(j-1)*m)=g3(y(j));  i=1;
  %with power input from left side
  i=1;
  A(i+(j-1)*m,i+(j-1)*m)=-3/(2*h);
  A(i+(j-1)*m,i+1+(j-1)*m)=4/(2*h);
  A(i+(j-1)*m,i+2+(j-1)*m)=-1/(2*h);
  b(i+(j-1)*m)=-P/(L*delta*K); 
end

for j=(n-1)/2+1:n-1;
    i=1;
    A(i+(j-1)*m,i+(j-1)*m)=-3/(2*h)+H/K;
    A(i+(j-1)*m,i+1+(j-1)*m)=4/(2*h);
    A(i+(j-1)*m,i+2+(j-1)*m)=-1/(2*h);
    b(i+(j-1)*m)=0;
end
%   i=2;
%   A(i+(j-1)*m,i+(j-1)*m)=4/(2*h);
%   b(i+(j-1)*m)=-P/(L*delta*K);
%   i=3;
%   A(i+(j-1)*m,i+(j-1)*m)=-1/(2*h);
%   b(i+(j-1)*m)=-P/(L*delta*K);
  
%   i=m-2;
%   A(i+(j-1)*m,i-2+(j-1)*m)=-1/(2*h);
%   b(i+(j-1)*m)=0; 
%   i=m-1;
%   A(i+(j-1)*m,i+(j-1)*m)=4/(2*h);
%   b(i+(j-1)*m)=0; 
for j=2:n-1;
  i=m;
  A(i+(j-1)*m,i+(j-1)*m)=-3/(2*h)+H/K;
  A(i+(j-1)*m,i-1+(j-1)*m)=4/(2*h);
  A(i+(j-1)*m,i-2+(j-1)*m)=-1/(2*h);
  b(i+(j-1)*m)=0; 
end 
v=A\b; % solve for solution in v labeling 
v=v+u_b;
w=reshape(v(1:mn),m,n); %translate from v to w 
mesh(x,y,w') 

