function [] = Laplace2d(N,M,L)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Laplace2d solves Laplace equation using FDM scheme. N: no of x axis mesh intervals, M: no of y axis mesh intervals
% It also compares the solutions for the equation:
% u_xx+u_yy=0, u(x,L)=2*sin(pi*x/L), u(x,0)=0=u(0,y)=u(L,y)

% Exact solution: u(x,y)=2*sin(pi*x/L)*sinh(pi*y/L)/sinh(pi)

% Run with typing in command window: a) Laplace2d(100,100,1) b) Laplace2d(50,50,4) c) Laplace2d(10,10,0.2)
% d) Laplace2d(1000,100,1)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

dx=L/N
dy=L/M
xx=0:dx:L;
yy=0:dy:L;

uu=@(x,y) 2*sin(pi*x/L).*sinh(pi*y/L)/sinh(pi);
[XX,YY]=meshgrid(xx,yy);
U=uu(XX,YY);

A=A2d(N,M,L);
%full(A)



u=zeros(M+1,N+1);
u(M+1,:)=2*sin(pi*xx/L);

temp=u;
temp([1 end],:) = [];
temp(:,[1 end]) = [];
temp(M-1,:)=temp(M-1,:)-u(M+1,2:N)/dy/dy;

T=A\temp(:);
u(2:M,2:N)=reshape(T,M-1,N-1);


mesh(xx,yy,u)

xlabel('x');
ylabel('y');
title('Solution of Laplace equation')
Error_Laplace =norm(u-U,Inf)

function A=A1d(a,b,J)
% A1D one dimensional finite difference approximation
% A=A1d(a,b,J) computes a sparse finite difference
% approximation of the one dimensional operator Delta on the
% domain Omega=(a,b) using J interior points
h=(b-a)/J;
e=ones(J-1,1);
A=spdiags([e/h^2 -(2/h^2)*e e/h^2],[-1 0 1],J-1,J-1);

function A=A2d(nx,ny,L)
% A2D finite difference approximation of -Delta in 2d
% A=A2d(nx,ny,L); constructs the finite difference approximation to
% -Delta on a nx x ny grid with spacing h=L/ny and homogeneous
% Dirichlet conditions all around
h=L/ny; Dxx=A1d(0,L,nx); Dyy=A1d(0,L,ny);
A=kron(speye(size(Dxx)),Dyy)+kron(Dxx,speye(size(Dyy)));