I don't get an error when running your code with MATLAB version 2019b, but the bigger problem you have is that your method will not scale well. A few tips off the top:
- Don't use anonymous functions unless absolutely necessary; they are not needed in this problem!
- Vectorize. You don't need to construct A & b in a for loop.
- Use spdiags. For a 5-point stencil there are 5 non-zero elements in each row (forming 5 diagonals), and you can use spdiags to generate your sparse matrix from an Nx5 matrix, rather than using spalloc and filling entries manually.
- Use exp(), not e^
I put together the following solution, which runs 
values 50, 100, and 500 in 0.027, 0.11, ands 1.56 seconds respectively on my computer. Your code runs the same values in 0.074, 0.31, and 405 seconds. That's >250 times faster for N=500!
values 50, 100, and 500 in 0.027, 0.11, ands 1.56 seconds respectively on my computer. Your code runs the same values in 0.074, 0.31, and 405 seconds. That's >250 times faster for N=500!close all
clear all
clc
% Grid Definition
nx = 50;
N = nx*nx;
x = linspace(0,1,nx);
y = x;
[X,Y] = ndgrid(x,y);
dx = x(2)-x(1);
% Compute RHS with no BC, assuming f(X,Y)
b = -1.36364*X.^4.5 + X.^3.5 - 0.818182*X.^9 + 0.818182*X.^8 - 0.0808081;
b = b * 247.5 .* (Y - 1).*Y.*cos(X);
b = b -45 *(-X.^5.5 + X.^4.5 - 0.444444*X + 0.222222).*(Y - 1).* Y.* sin(X);
b = b -20 *(X - 1).*X.*cos(X) ;
b = exp( X.^4.5 ) .* b;
b = b.*dx.*dx;
% Overwrite BC's on the RHS using vectorized assignment
isDirichlet = (X==0) | (X==1) | (Y==0) | (Y==1);
b(isDirichlet) = 0.0;
% Compute A Matrix diagonals Ad(1:N,1:5) assuming normal stencils
Ad = zeros(N,5);
% Compute off-center points, note that there is no y dependence in
% k(x,y)=cos(x)
Ad(:,1) = -cos(X(:)); % Down = -kd = -cos(X)
Ad(:,2) = -cos(X(:)-dx*0.5); % Left = -kl = -cos(X-dx/2)
Ad(:,4) = -cos(X(:)+dx*0.5); % Right= -kl = -cos(X+dx/2)
Ad(:,5) = Ad(:,1); % Up = -ku = -cos(X) = Down
% Use these to get center by summing the negatives
Ad(:,3) = -sum( Ad(:,[1,2,4,5]), 2 ); % Center
% Overwrite BC's on the A Matrix diagonals
idx = find(isDirichlet(:)); % array of inidices of BC points
Ad(idx,:) = 0.0; % Zero all elements of rows defined by idx
Ad(idx,3) = 1.0; % Make center node=1
% Create A Matrix from Ad, assigning columns to diagonals
A = spdiags(Ad,[-nx,-1,0,1,nx],N,N)'; % ADDED TRANSPOSE
% Solve & Reshape
uh=A\b(:);
uh = reshape(uh,[nx,nx]);
% Compute Exact solution and error
u_exact = 10*X.*Y.*(1-X).*(1-Y).*exp(X.^(4.5));
u_error = abs( (uh - u_exact)./u_exact );
% Reshape and plot solution on two subplots
subplot(1,3,1); view(0,90); hold on; colorbar;
surf(x,y,uh,'edgecolor','none');
xlabel('x'); ylabel('y')
title('Computed')
shading interp
axis square
% Add exact solution
subplot(1,3,2); view(0,90); hold on; colorbar;
surf(x,y,u_exact,'edgecolor','none');
xlabel('x');
title('Exact')
shading interp
axis square
% add Error
subplot(1,3,3); view(0,90); hold on; colorbar;
surf(x,y,u_error,'edgecolor','none');
xlabel('x');
title('Error')
shading interp
axis square
