function [A_BPA, b_BPA] = calc_stiffness_matrix_BPA( grid, ice)
% Calculate stiffness matrix and load vector for the BPA

% Assemble the stiffness matrix for Eq. 1

% d/dx ( 2 eta ( 2 du/dx + dv/dy ))
Au1 = grid.M_ddx_ak_bk  * 2 * diag_sparse( ice.eta_ak ) * (2 * grid.M_ddx_bk_ak  * grid.M_map_ijku2ijk + grid.M_ddy_bk_ak * grid.M_map_ijkv2ijk);
% d/dy (   eta (   du/dy + dv/dx ))
Au2 = grid.M_ddy_ak_bk  *     diag_sparse( ice.eta_ak ) * (    grid.M_ddy_bk_ak  * grid.M_map_ijku2ijk + grid.M_ddx_bk_ak * grid.M_map_ijkv2ijk);
% d/dz (   eta (   du/dz         ))
Au3 = grid.M_ddz_bks_bk *     diag_sparse( ice.eta_bks) *      grid.M_ddz_bk_bks * grid.M_map_ijku2ijk;
% Combine
Au = grid.M_map_ijk2ijku * (Au1 + Au2 + Au3);

% Assemble the stiffness matrix for Eq. 2

% d/dy ( 2 eta ( 2 dv/dy + du/dx ))
Av1 = grid.M_ddy_ak_bk  * 2 * diag_sparse( ice.eta_ak ) * (2 * grid.M_ddy_bk_ak  * grid.M_map_ijkv2ijk + grid.M_ddx_bk_ak * grid.M_map_ijku2ijk);
% d/dx (   eta (   dv/dx + du/dy ))
Av2 = grid.M_ddx_ak_bk  *     diag_sparse( ice.eta_ak ) * (    grid.M_ddx_bk_ak  * grid.M_map_ijkv2ijk + grid.M_ddy_bk_ak * grid.M_map_ijku2ijk);
% d/dx (   eta (   dv/dz         ))
Av3 = grid.M_ddz_bks_bk *     diag_sparse( ice.eta_bks) *      grid.M_ddz_bk_bks * grid.M_map_ijkv2ijk;
% Combine
Av = grid.M_map_ijk2ijkv * (Av1 + Av2 + Av3);

% Combine Eqs. 1 and 2
A_BPA = Au + Av;

% Fill load vector
b_BPA = zeros( grid.nijkuv, 1);
for i = 1: grid.nx
for j = 1: grid.ny
for k = 1: grid.nz
  rcu = grid.ijkuv2rc( i,j,k,1);
  rcv = grid.ijkuv2rc( i,j,k,2);
  b_BPA( rcu) = ice.rho * ice.g * ice.dh_dx_b( i,j);
  b_BPA( rcv) = ice.rho * ice.g * ice.dh_dy_b( i,j);
end
end
end

end