Use none of them. This is NOT an ODE problem, except that you seemed to want to write it as such. It is certainly not something that dsolve will do, because there is nothing symbolic about the problem.
All you are looking to do is integrate a spline.
fnint will do that, as long as you have the curve fitting toolbox. You want to integate twice.
Because you want to set up a boundary value at x == 3, we need to recognize what is the integral of a constant, and use that to determine the "initial condition" to use for the inner integral. Remember that the value for the second argument to fnint is just a constant thjat will be added to the entire spline, a constant offset. Think of it as the constant of integration.
>> MIItemp = fnint(fnint(M));
>> y_3 = fnval(MIItemp,3)
>> MII = fnint(fnint(M,-y_3/3),0);
MII is itself a spline.
struct with fields:
breaks: [1×31 double]
coefs: [30×6 double]
You can evaluate it, plot it, whatever. In fact, this is the analytical indefinite double integral of the function M, subject to the "boundary" values you indicated. It is defined as a spline, because the original function M was also defined as a spline.