2D wave equation with Mur boundary condition - setting up the RHS and solving (time-steps)

Question

I am trying to solve a 2D wave equation implicitly using FD with central approximations with the following boundary conditions

$$\begin{align} &u=2\sin\left(\frac{2\pi}{5}t\right)\quad \text{at }x=0\\ &\frac{\partial u}{\partial y}=0 \quad \text{at }y=0\\ &\frac{\partial u}{\partial y}=0 \quad \text{at }y=2\\ &\frac{\partial u}{\partial t}=-\frac{\partial u}{\partial x} \quad \text{at }x=5 \end{align}$$

graphically represented as

The boundary condition at x=0 generates the wave. The boundary condition at x=5 refers to Mur boundary condition, i.e. the wave doesn't bounce back, but simply continues to move outside the domain. The general equation is given as

The Mur boundary condition can be expressed mathematically as

The boundary conditions were evaluated analytically and these are the 7 equations following the entire domain:

Problem- I have managed to write some code, including setting up the [A] matrix (the LHS of the equations), but I have problems with formulating the [b] matrix (RHS) and solving it through the time-steps. I would appreciate if someone show me how to write it for the first equation - then I will be able to do it dor the remaining 6.

% length, time, height
L = 5; % [m]
h = 2; % [m]

d = 0.25; % space spacing
dt = 0.05; % time increment
M = 100; 
T_max = M*dt; % [s]

nx = floor(L/d); % number of <x> samples
ny = floor(h/d); % number of <y> samples
k = floor(T_max/dt) + 1; % number of time samples

% Constants
alpha = dt/d;
r = dt^2/d^2;

% Number of grid points:
N=nx*ny;

% Initialise matrices
U = zeros(N,1);
b = zeros(N,1);
x = zeros(N,1);
A = zeros(N,N);

% Set numbering for the [A] matrix
num = 1;
for i=1:nx
    for j=1:ny
        number(i,j) = num;
        num = num + 1;
    end 
end 

% [A] matrix - boundary conditions
    % Case 1 
    for i=2:nx-1
        for j=2:ny-1
            ii=number(i,j);
            A(ii, number(i, j)) = 1+4*r;
            A(ii, number(i+1, j)) = -r;
            A(ii, number(i-1,j)) = -r;
            A(ii, number(i, j+1)) = -r;
            A(ii, number(i, j-1)) = -r;
        end 
    end       

    % Case 2
    for i=1
        for j=2:ny-1
            ii=number(i,j);
            A(ii, number(1, j)) = 1+4*r;
            A(ii, number(2, j)) = -r;
            A(ii, number(1, j+1)) = -r;
            A(ii, number(1, j-1)) = -r;
        end 
    end 

    % Case 3 
    for i=nx
        for j=1:ny
            ii=number(i,j);
            A(ii, number(nx, j)) = 1+alpha;
            A(ii, number(nx-1,j)) = 1-alpha;
        end 
    end 

    % Case 4 
    for i=2:nx-1
        for j=1
            ii=number(i,j);
            A(ii, number(i, j)) = 1+4*r;
            A(ii, number(i-1,j)) = -r;
            A(ii, number(i+1,j)) = -r;
            A(ii, number(i, j+1)) = -2*r;
        end 
    end 

    % Case 5
    for i=1
        for j=1
            ii=number(i,j);
            A(ii, number(i,j)) = 1+4*r;
            A(ii, number(i+1,j)) = -r;
            A(ii, number(i, j+1)) = -2*r;
        end 
    end 

    % Case 6
    for i=2:nx-1
        for j=ny
            ii=number(i,j);
            A(ii, number(i,j)) = 1+4*r;
            A(ii, number(i-1,j)) = -r;
            A(ii, number(i+1, j)) = -r;
            A(ii, number(i, j-1)) = -2*r;
        end 
    end 

    % Case 7
    for i=1
        for j=ny
            ii=number(i,j);
            A(ii, number(i,j)) = 1+4*r;
            A(ii, number(i+1, j)) = -r;
            A(ii, number(i, j-1)) = -2*r;
        end
    end

% Inverse of [A] matrix
A_inv = inv(A);

% [b] matrix (RHS)
time = 0;
for t=1:M
    % Computing [b] matrix
    for i=1:nx
        for j=1:ny
            ii = number(i,j);
            % U(1,1) = 2*sin((2*pi/5)*(t*dt));

            % Case 1
            for i=2:nx-1
                for j=2:ny-1

                end 
            end 
            % or

            for i=1:nx
                for j=1:ny

                end 
            end 

        end 
    end 

    % The solution vector
    x = A_inv*b;

    % 2D matrix from the solution vector