clear all; close all

global m h uxc

m=20; h=1/m; x = 0:h:1; uxc = 5; 
y=x; u0=zeros(m+1,m+1); f=u0; ue=u0; h1 = h*h;

% Set the right hand side

for i=1:m+1,
  for j=1:m+1
    f(i,j) = ((1-pi^2) + uxc)*bc(x(i),y(j));
    ue(i,j) = bc(x(i),y(j));
  end
end

% Adjust the right hand side

for i=2:m,
   f(i,2) = f(i,2) - bc(x(i),y(1))/h1;
   f(i,m) = f(i,m) - bc(x(i),y(m+1))/h1;
end

for j=2:m,
  f(2,j) = f(2,j) - bc(x(1),y(j))/h1 + uxc*bc(x(1),y(j))/(2*h);
  f(m,j) = f(m,j) - bc(x(m+1),y(j))/h1 - uxc*bc(x(m+1),y(j))/(2*h);
end

rhs = u2d_to_1d(m,f);

mt = length(rhs); x0 = zeros(mt,1); 

params=[1.0e-8,1000,1];

[xs, error, total_iters] = my_gmres(x0, rhs, 'atv', params);

u = u1d_to_2d(m,xs);

for i=1:m+1,                    % Add boundary conditions
  u(i,1) = bc(x(i),y(1));
  u(i,m+1) = bc(x(i),y(m+1));
end

for j=2:m,
  u(1,j) = bc(x(1),y(j));
  u(m+1,j) = bc(x(m+1),y(j));
end
  
figure(1); mesh(x,y,u);xlabel('x');ylabel('y');
zlabel('u'); title('Solution Plot');
figure(2); mesh(x,y,u-ue); xlabel('x');ylabel('y');
zlabel('e(x,y)'); title('Error Plot');
total_iters,  norm(u-ue,'inf')