lptest.m
format long
err_ul = [];
err_ur = [];
err_c = [];
err_ll = [];
err_lr = [];
for ind=3:6

appr = poisson13(0,1,0,1,2^ind,2^ind);

clos = infser(2^ind,2^ind);

diff = appr(2^ind+1,2^ind+1) - clos(2^ind+1,2^ind+1);

err_ul = [err_ul; appr(1,1) - clos(1,1)];
err_ur = [err_ur; appr(1,2^ind+1) - clos(1,2^ind+1)];
err_c = [err_c; appr(2^(ind-1)+1,2^(ind-1)+1) - clos(2^(ind-1)+1,2^(ind-1)+1)];
err_ll = [err_ll; appr(2^ind+1,1) - clos(2^ind+1,1)];
err_lr = [err_lr; appr(2^ind+1,2^ind+1) - clos(2^ind+1,2^ind+1)];

%err = [err; max(max(diff))];
end

err_ul 
err_ur
err_c
err_ll
err_lr

poisson13.m

% Program 8.5 Finite difference solver for 2D Poisson equation 
%   with Dirichlet boundary conditions on a rectangle 
% Input: rectangle domain [xl,xr]x[yb,yt] with MxN space steps 
% Output: matrix w holding solution values 
% Example usage: w=poisson13(0,1,0,1,10,10) 
function w=poisson13(xl,xr,yb,yt,M,N) 
T0 = 0;
T1 = 10;
h=(xr-xl)/M;h2=h^2;k=(yt-yb)/N;k2=k^2;

f=@(x,y) 0; %; % define input function data 
g1=@(x) 0; % define boundary values 
g2=@(x) 0; % Example 8.8 is shown 
g3=@(y) 0;
g4=@(y) 10; 
m=M+1;n=N+1; mn=m*n; 
 
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; 
    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)=f(x(i),y(j)); 
  end 
end 
for i=1:m % bottom and top boundary points 
  j=1;A(i+(j-1)*m,i+(j-1)*m)=1;b(i+(j-1)*m)=g1(x(i)); 
  j=n;A(i+(j-1)*m,i+(j-1)*m)=1;b(i+(j-1)*m)=g2(x(i)); 
end 
for j=2:n-1 % left and right boundary points 
  i=1;A(i+(j-1)*m,i+(j-1)*m)=1;b(i+(j-1)*m)=g3(y(j)); 
  i=m;A(i+(j-1)*m,i+(j-1)*m)=1;b(i+(j-1)*m)=g4(y(j)); 
end 
v=A\b; % solve for solution in v labeling 
w=reshape(v(1:mn),m,n); %translate from v to w 

infser.m

function A=infser(N,M)
A=rv(N, M, 100)';

function msol=rv(N, M, st)
h = 1/N;
k = 1/M;

msol = zeros(M+1,N+1);
for i=1:(M+1)
    for j =1:(N+1)    
        msol(i,j) = infsum(st, 0+ h*(i-1), 0+ k*(j-1));
    end
end

function cv=infsum(iter, x, y)

cv = 0;

for k=0:(iter-1)
    cv = cv + (sinh((2*k+1)*pi*y)/sinh((2*k+1)*pi)) * (40/((2*k+1)*pi))* sin((2*k+1)*pi*x);
    
end