The problem is also present in the development version; we must care to not forget the mass matrix when computing the solution. In ode45.m, line 262 :

 if (! isempty (odeopts.Mass) && isnumeric (odeopts.Mass))
    havemasshandle = false;
    mass = odeopts.Mass;  # constant mass
  elseif (isa (odeopts.Mass, "function_handle"))
    havemasshandle = true;    # mass defined by a function handle
  else  # no mass matrix - creating a diag-matrix of ones for mass
    havemasshandle = false;   # mass = diag (ones (length (init), 1), 0);
  endif

the mass matrix is defined in the constant case but not used. I suggest to use the parameter "fun" the same way a dynamic mass matrix is handled, e.g. line 301, add the bottom lines :

if (havemasshandle)   # Handle only the dynamic mass matrix,
    if (massdependence) # constant mass matrices have already
      mass = @(t,x) odeopts.Mass (t, x, odeopts.funarguments{:});
      fun = @(t,x) mass (t, x, odeopts.funarguments{:}) ...
             \ fun (t, x, odeopts.funarguments{:});
    else                 # if (massdependence == false)
      mass = @(t) odeopts.Mass (t, odeopts.funarguments{:});
      fun = @(t,x) mass (t, odeopts.funarguments{:}) ...
             \ fun (t, x, odeopts.funarguments{:});
    endif

# here the suggested fix

else
     invmass = inv(mass);
     fun = @(t,x) invmass * fun (t, x, odeopts.funarguments{:});
endif

Do not forget to uncomment the mass matrix line 268 (it will be equal to identity)

Package : odepkg v0.8.5
Functions : any solver (ode45, ode54, etc) used in the package

Description :
Any use of the 'Mass' option with one of the solver (though the settings function 'odeset') will fail if the mass matrix is constant; the function runs without error/warning but the result is incorrect.

This problem does not occur when a function handle @mass is passed as the mass option.

Fix :
The change must be done in every solver. For example in ode54.m, line 232, mass definition must be uncommented :

vhavemasshandle = false; %# vmass = diag (ones (length (vinit), 1), 0);

to become :

vhavemasshandle = false;
vmass = diag (ones (length (vinit), 1), 0);

and then used line  367 :

vk(:,j) = feval ...
          (vfun, vthetime, vtheinput, vfunarguments{:});

must be changed into

vk(:,j) = vmass \ feval ...
          (vfun, vthetime, vtheinput, vfunarguments{:});

Important remark : in the case of a constant matrix we dont need to perform a matrix division at each solver step, because the computation time will greatly increase. It is better to compute the inverse of the mass matrix before the main loop, and then do a multiplication (rather than a matrix division). This gives :

l.227 :

if (~isempty (vodeoptions.Mass) && isnumeric (vodeoptions.Mass))
    vhavemasshandle = false;
    invmass = inv(vodeoptions.Mass); %# constant mass
  elseif (isa (vodeoptions.Mass, 'function_handle'))
    vhavemasshandle = true; %# mass defined by a function handle
  else %# no mass matrix - creating a diag-matrix of ones for mass
    vhavemasshandle = false;
    invmass = diag (ones (length (vinit), 1), 0);
  end

and line 356 :

if (vhavemasshandle)
...
...
...
        vk(:,j) = vmass \ feval ...
          (vfun, vthetime, vtheinput, vfunarguments{:});
      else
        vk(:,j) = invmass * feval ...
          (vfun, vthetime, vtheinput, vfunarguments{:});
      end