lptest.m
N = 64;
M = 64;
h = 1/N;
k = 1/M;
x = (0:M)*h;
y = (0:N)*k;

A=rv(N, M, 100);
mesh(x',y',A)

rv.m

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

infsum.m

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