Integral for Biot-Savart Law on a box

Question

I'm trying to make a basic computer model of a bar magnet. In the process I came across this question and answer that appears to have an appropriate equation for me to use. I say appears because my maths knowledge is only barely past high-school level, and the answerer stops here:

you can simply use the Biot-Savart law to calculate the magnetic field: $$\mathbf B(\mathbf x) = \frac{\mu_0}{4\pi}\int_{\mathbb S}d\mathbf a' \ \mathbf K(\mathbf x') \times \frac{\mathbf{x-x'}}{|\mathbf{x-x'}|^3}$$ I believe you can take it from here.

So I understand that in this case I can use the fact that the net magnetic field (the thing I want to model) is composed of the sum of all the magnetic fields produced in this situation. Which is to say I need to add up the field produced by each of the faces with respect to each point I want to model.

I also figure that that is what the part of the equation that I don't understand is trying to express, namely the section $\int_{\mathbb S}d\mathbf a' \ \mathbf K(\mathbf x')$

So in the above example a bar magnet is being modeled like so:

You can model a bar magnet by a rectangular box with a constant magnetization in one direction. Let's take the box $[0,a]\times[0,b]\times[0,c]$, with a constant magnetization $\mathbf M(\mathbf x) = M_0 \ \hat{\mathbf k}$, where $\hat{\mathbf k}$ is the unit vector in the $z$ direction. The bound volume and surface current densities are: $$\mathbf J_b(\mathbf x) = \boldsymbol{\nabla}\times\mathbf M(\mathbf x)$$ $$\mathbf K_b(\mathbf x) = \mathbf M(\mathbf x) \times \hat {\mathbf n}$$ The volume current density is zero because $\mathbf M$ is constant. For the surface current density, the top and bottom faces don't contribute since $M_0 \hat{\mathbf k}\times\hat {\mathbf k}=0$. For the other four faces we have: $$\mathrm{x=0 \ face:} \ \mathbf K_1 = M_0 \ \hat{\mathbf k}\times (-\hat{\mathbf i}) = -M_0 \ \hat{\mathbf j}$$ $$\mathrm{x=a \ face:} \ \mathbf K_2 = M_0 \ \hat{\mathbf k}\times \hat{\mathbf i} = M_0 \ \hat{\mathbf j}$$ $$\mathrm{y=0 \ face:} \ \mathbf K_3 = M_0 \ \hat{\mathbf k}\times (-\hat{\mathbf j}) = M_0 \ \hat{\mathbf i}$$ $$\mathrm{y=b \ face:} \ \mathbf K_4 = M_0 \ \hat{\mathbf k}\times \hat{\mathbf j} = -M_0 \ \hat{\mathbf i}$$ Now that you know the bound current distribution, you can simply use the Biot-Savart law to calculate the magnetic field: $$\mathbf B(\mathbf x) = \frac{\mu_0}{4\pi}\int_{\mathbb S}d\mathbf a' \ \mathbf K(\mathbf x') \times \frac{\mathbf{x-x'}}{|\mathbf{x-x'}|^3}$$ I believe you can take it from here.

My question is how do I evaluate the integral portion of this equation? I'm looking to turn this into a piece of computer code, and my background is pretty shallow when it comes to this level of maths.

Edit: I understand Matlab has an integrate function, but I would prefer not to buy a license for that if possible.

Edit2:

After thinking about this some more, and with the help of Ian's comments I have determined what I think I need to do, which is best expressed graphically by the diagram I have just drawn:

Answer 1

So it seams to me the specific translation you were seeking looks something like this, and probably would have better been asked on physics or stack overflow.

for each position (x) in space:

    for each face in the set of faces:
        
        #create a function for this faces surface current
        Kface(dx) = M_0 cross surface_normal ??? This function feels more like a simple vector
    
        for each point (df) on current face:
        
             B_x_df = mu_0 / 4 * pi * Kface(dx) cross x-dx / mag(x-dx)^3

    B_x = sum of intermediate B_x_df values
```