



function [U,XX,YY]=fdmdir2D16(f,a,b,N,c)
%fdm solver for -Laplacian u+cu=f in (a,b)x(a,b) (c>=0 nonnegative constant)
%u(s)=0 on bnd
%INPUT
%f - function handle to rhs f (it must work for vectors i.e .f(x)=[f(x(1),...,f(x(N-1))]
%a,b- interval ends
%c - positive constant parameter of PDE
%N- no od mesh points minus one in each direction
%OUTPUT
%U - solution as a matrix i.e. U(k,l)=u(l*h,k*h)
%XX,YY - meshgrid marices (XX constant cols, YY constant rows)
%EXAMPLES
%s=@(x,y) sin(x).*sin(y);f=@(x,y)2*s(x,y); c=0; a=0;b=pi;  N=20; [U,XX,YY]=fdmdir2D16(f,a,b,N,c);mesh(ER=s(XX,YY)-U);er0=norm(reshape(ER,(N-1)*(N-1),1),'inf');pause(1);N=2*N;[U,XX,YY]=fdmdir2D16(f,a,b,N,c);mesh(ER=s(XX,YY)-U);[er0,er1=norm(reshape(ER,(N-1)*(N-1),1),'inf'),er0/er1]

h=(b-a)/N;
n=(N-1)*(N-1);
%3 diagonal matrix N-1 X N-1 with 2 on main diag and -1 on sub super diags
L1=2*sparse(diag(ones(N-1,1))) - sparse(diag(ones(N-2,1),1)) - sparse(diag(ones(N-2,1),-1));
I1D=speye(N-1);
A=kron(I1D,L1)+kron(L1,I1D)+ h*h*c*speye(n);

%mesh
x=linspace(a,b,N+1); x=x';

%h scaling goes also into rhs .e. F
[XX, YY] = meshgrid(x(2:N));
FF=h*h*f(XX,YY);
F=reshape(FF,n,1);
u=A\F;
U=reshape(u,N-1,N-1);

end

