second sensitivity equations for large scale ODEs with multiple parameters

In order to calculate geodesic path,  I will calculate  second sensitivity equations for 19 ODEs with 41 parameters. Owing to integrate for the second  sensitivity equations, the equations should be keep and I use symbolic expression and matlabfunction to do it. But the matlab shows me the error information: out of memory for the second step( as shown next)

1. [SS,dSS,SS0] = sensitivity_equations_another([obj.x;obj.S(:)], obj.theta, [obj.dx;obj.dS(:)], [obj.x0;obj.S0(:)],M);

2 obj.ODE_2nd_order_sensitivity_rhs=matlabFunction([obj.dx;obj.dS(:);obj.dSS(:)],'vars', {'t',[obj.x;obj.S(:);obj.SS(:)], obj.theta});

3. ode45(@(t,x_values)obj.ODE_2nd_order_sensitivity_rhs(t,x_values,para),obj.time_points_to_evaluate,obj.ODE_2nd_order_sensitivity_init(para));

here the size of obj.x is 19*1; the size of obj.S() is 798*1 and the sizie of obj.SS(:) is 33516*1

even if the matrix was seperated for different parts, the error information still is : the symbolic expression is too large to calculate

Then for the large scale ODEs and multiple parameters, how to solve it