function error = TwoD_Case(N,omega)
    h = 1/N;
    nu = 0.1;
    x = [0:h:1]';
    y = [0:h:1]';
    % Initial value and boundary condition
    u = zeros(N+1);
    ue = zeros(N+1);
    for j = 2:N
        for i = 2:N
            u(i,j) = sin(pi*x(i))*sin(pi*y(j));
        end
    end
    % Solve the linear system iteratively by SOR method
    kmax = 1000; 
    tol = 1.0E-13; 
    f = zeros(N+1);
    % Solve for u from t = 1 to t = N
    for t = 1:N
        for j = 2:N
            for i = 2:N
                f(i,j) = (nu/2)*(u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1))+(h-2*nu)*u(i,j);
            end
        end
        k = 0;
        rmax = 1;
        while (k<kmax)&&(rmax>tol)
            rtmp = 0;
            for j = 2:N
                for i = 2:N
                    % The initial u will be the prior u we have solved
                    uGS = ((nu/2)*(u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1))+f(i,j))/(h+2*nu);
                    utmp = u(i,j)+omega*(uGS-u(i,j));
                    rtmp = max(rtmp, abs((h+2*nu)*(utmp-u(i,j))));
                    u(i,j) = utmp;
                end
            end
            rmax = rtmp;
            k = k+1;
        end
    end
    % The exact solution
    for j = 2:N
        for i = 2:N
            ue(i,j) = exp(-2*nu*pi^2)*sin(pi*x(i))*sin(pi*y(j));
        end
    end
    % The error of approximation
    error = 0;
    for j = 2:N
        for i = 2:N
            error = max(error,abs(ue(i,j)-u(i,j)));
        end
    end
end