Solving 2D Numerical Quadrature

Right now, I am trying to code for 2D Left-Point Rectangular Rule for this equation:

The equation that I am supposed to use is provided (I^LP)

where I_i,j is a double summation (the furtherest right equation)

I know now that i can't increase both variables at the same time, but I can't figure how to do it properly even after I tried nested for loops. I am too burnt out to think properly anymore so if someone could help me, it would be greatly appreciated.

[J. Alex Lee]

show your attempted code for nested loop - that would be the most straightforward way

[Dominic Marquez]

This is what I have for 2D left point rectangular rule, the problem I am having is that my answers for 5, 10, and 20 subintervals are so far off the analytical answer of 21.333 that makes me think I'm doing it wrong. However, if you use 300 + subintervals you end up with a close answer, so I don't know if this is a quirk of rectangular rule or if I'm doing it wrong.

a = 0; b = 4; 
c = -2; d = 2;
sum = 0;

prompt1 = ['Please choose an integration scheme: \n(1) Left-Point Rectangular Rule \n' ...
    '(2) Right-Point Rectangular Rule \n(3) Mid-Point Rectangular Rule \n' ...
    '(4) Trapezoidal Rule \n(5) Simpson''' 's 1/3 Method'];
selection = input(prompt1);

if selection < 1 || selection > 5
    error('Please choose a valid integration scheme.')
end

prompt2 = 'How many subintervals are there in the x-direction?';
n_x = input(prompt2);
prompt3 = 'How many subintervals are there in the y-direction?';
n_y = input(prompt3);

h_x = (b - a) / n_x;
h_y = (d - c) / n_y;

if selection == 1
    for y = c : h_y : (d - h_y)
        for x = a : h_x : (b - h_x)
            sum = sum + (x^2 - 3 * y^2 + x * y^3);     
        end
    end
    
lp = sum * h_x * h_y;
disp(['Left-Point Rectangular Rule = ' num2str(lp)]);

I have the other integration schemes if needed.

[J. Alex Lee]

How far off is "so far off", and what reason do you have to expect closer answers? The implementation looks right.

Is the purpose of this exercise to understand the numerical quadrature methods? Or how to implement the methods in code? If it's to understand the method, one exercise you can do is plot the error from true solution as a function of n_x (for simplicity set n_y=n_x), and notice the relationship. Then do the same thing for trapezoidal rule and simpson's rule.