function test_ode15s_jac m = 3.0; k = 1.0; naturalOmega = sqrt(k/m); massMat=[1 0; 0 m]; options=odeset('Stats','on'); tol = 1e-6; options=odeset(options, 'RelTol', tol); options=odeset(options, 'AbsTol', tol); jacFunc = @(time, u) calcJacobianMat(time, u, k); options=odeset(options, 'Jacobian', jacFunc); options=odeset(options, 'Mass', massMat); period = 2*pi/naturalOmega; n=3; tlist = linspace (0, n*period, n*200); y0=[1 0]'; yAnal=analSoln(tlist, m, k); rhsFunc = @(t, u) calcRHSODE(t,u, k); [outTimes,u2]=ode15s(rhsFunc, tlist, y0, options); maxErr=max(abs(u2(:)-yAnal(:))); fprintf('tol=%12.3e, maxErr=%12.3e\n', tol, maxErr); if(0) figure; plot(outTimes, u2(:,1)); %plot(u2(:,1), u2(:,2)); grid; end end function rhsODE=calcRHSODE(t, u, k) jac=[0 -1; k 0]; rhsODE=-jac*u; end function jac=calcJacobianMat(time, u, k) jac=[0 1; -k 0]; end function mm=calcMassMat(time, u, m) mm=[1 0; 0 m]; end function y=analSoln(t, m, k) omega = sqrt(k/m); v = -omega*sin(omega*t); u = cos(omega*t); y=[u;v]'; end