Solve Laplace equation in 2D with nonhomogeneous boundary conditions, using Galerkin method

Question

I want to solve Laplace equation:

$$\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}=0$$

in the domain $0<x<1$, $0<y<1$, with non-homogeneous boundary conditions $u(x=0,y)=1$ and $u(x=1,y)=u(x,y=0)=u(x,y=1)=0$.

An analytical solution using Fourier series is simple and gives the solution:

Now I want to (approximately) solve it using Galerkin method.

I am following Fletcher's book (Fletcher - Computational Galerkin Methods). I have checked here and here.

I tried writing a code to find the approximate solution, given in the form

$$F_{ij}(x,y)=\sin(i \pi x) \sin(j \pi y)$$

(I tried also with cos, and sin+cos).

and the initial approximation $u_0(x,y) = u_0 = (1-x)$ to satisfy the B.C. Solution should be

$$ u\approx \sum_i \sum_j a_{ij} F_{ij} + u_0$$

The problem: since the basis function are orthogonal to each other, by the method of weighted residuals, all coefficients are always zero.

I even tried to compute it with the code:

clc; clear;

syms x y real positive
syms i j integer positive

% differential equation
L = @(u) diff(u,x,2) + diff(u,y,2); % linear operator on differential equation
% L = @(u) diff(u,x,2); % linear operator on differential equation
% f = sym(1); % R.H.S.
f = sym(0); % R.H.S.

% boundary location
xl = sym([0 1]);
yl = sym([0 1]);

pi = sym(pi)
F = sin(i*pi*x)*sin(j*pi*y) % trial function


u0 = (1-x)
M = 3;
N = 3;

ii = 1:M;
jj = 0:N;
ij = combvec(ii,jj);

for n=1:size(ij,2)
    phi(n) = subs(F,[i,j],[ij(1,n),ij(2,n)]); % basis functions
end

syms a0

a = sym("a",size(phi));
u = sym(0);
u = u+u0; % initial approximation, to satisfy BC
for n=1:length(a)
    u = u+a(n)*phi(n);
end
u;

% residual
r = L(u)-f+L(u0);

% inner product
inner = @(u,v) int(int(u*v,x,xl(1),xl(2)),y,yl(1),yl(2));

% form linear system of equations for coefficients aij
for n=1:length(a)
    R(n) = inner(r,phi(n)) == 0;
end
R'

a_coef = solve(R,a)

a_coef = struct2cell(a_coef)'
u_a = subs(u,a,a_coef);

I can't find what exactly is wrong here. Perhaps the trial functions are wrong?

How to solve this particular equation, with the given boundary conditions?