Jacobian of Massive Equations

As the title states I am trying to find the Jacobian of a few massive equations.

Per the documentation I have set up the problem as follows:

% Define symbols to be used throughout.
syms lat lon u v w h phi theta psi p q r delta1 delta2 delta3 delta4 Hg GM J2 Re e Ixx Ixz Iyy Izz omega m real
% Rotational rate in body axis.
omegabnb = [p;q;r];
% Inertia matrix.
I = [Ixx,0,Ixz;0,Iyy,0;-Ixz,0,Izz];
% Rotational rate of Earth.
omegaeie = [0;0;omega];
% Equatorial radius of curvature.
RE = Re/sqrt(1-e^2*sin(lat)^2);
% Polar radius of curvature.
RN = (Re*(1-e^2))/((1-e^2*sin(lat)^2)^(3/2));
% Euler orientations.
roll = [1,0,0;0,cos(phi),-sin(phi);0,sin(phi),cos(phi)];
pitch = [cos(theta),0,sin(theta);0,1,0;-sin(theta),0,cos(theta)];
yaw = [cos(psi),-sin(psi),0;sin(psi),cos(psi),0;0,0,1];
% Euler Matrix.
Cbn = simplify(yaw*pitch*roll);
Cnb = Cbn';
% LLA Matrix.
[~,Cne] = nRe(lat,lon);
Cne = simplify(Cne);
Cen = Cne';
% Earth rotation orientation.
[~,Cei] = eRi(Hg);
Cei = simplify(Cei);
Cie = Cei';
% body velocity
vbeb = [u;v;w];
vneb = Cbn'*vbeb;
vneb = simplify(vneb);
% Latitude and longitudinal rates.
latdot = simplify(vneb(1)/(RN+h));
longdot = simplify((vneb(2)/(RE+h))/cos(lat));
omeganen = [longdot*cos(lat);-latdot;-longdot*sin(lat)];
% Forces (outputs six equations determining the various forces and moments, equations in terms of most of the state)
% VERY long equations (126 coefficients).
[Fb,Mb] = MomentFunc;
% Radii calculations.
Rc = [1/(RN+h),0,0;
    0,1/((RN+h)*cos(lat)),0 ;
    0,0,-1];
% Position calculation.
reeb = [((Re/sqrt(1-e^2*(sin(lat)^2)))+h)*cos(lat)*cos(lon);
    ((Re/sqrt(1-e^2*(sin(lat)^2)))+h)*cos(lat)*sin(lon);
    ((Re*(1-e^2)/sqrt(1-e^2*(sin(lat)^2)))+h)*sin(lat)];
riib = Cie*reeb;
Riib = sqrt(riib(1)^2+riib(2)^2+riib(3)^2);
tmpx = (1-5*(riib(3)/Riib)^2)*riib(1);
tmpy = (1-5*(riib(3)/Riib)^2)*riib(2);
tmpz = (3-5*(riib(3)/Riib)^2)*riib(3);
tmp = [tmpx;tmpy;tmpz];
tmp = riib + (3/2)*J2*(Re^2)*tmp/(Riib^2);
% Acceleration due to gravity.
gammaiib = -GM/(Riib^3)*(tmp);
% Latitude, Longitude and Altitude rate.
gammadot = Rc*vneb;
% Body velocity rate.
vbebdot = Cen*Cnb*(Fb/m) + Cei*gammaiib - Cbn*Cne*SkewMat(Cen*omeganen+Cne*Cnb*omegabnb)*Cen*Cnb*vbeb;
% Angular rotational rate.
omegabnbdot = I\(Mb-SkewMat(Cbn*Cne*omegaeie+Cbn*omeganen+omegabnb)*I*(Cbn*Cne*omegaeie+Cbn*omeganen+omegabnb))-SkewMat(-omegabnb)*Cbn*omeganen;
% Roll rate.
phidot = p*q*sin(phi)*tan(theta)+r*cos(phi)*tan(theta);
% Pitch rate.
thetadot = q*cos(phi)-r*sin(phi);
% Yaw rate.
psidot = q*sin(phi)*sec(theta)+r*cos(phi)*sec(theta);
% Control deflection rate (DUE TO STATE, NOT DUE TO CONTROL).
deltadot = 0;
% Jacobian for linearization.
A = jacobian([gammadot;vbebdot;omegabnbdot;phidot;thetadot;psidot;deltadot;deltadot;deltadot;deltadot],[lat,lon,h,u,v,w,p,q,r,phi,theta,psi,delta1,delta2,delta3,delta4]);

This progam runs for several days and then ends with an out of memory error. I wish I'd stored the error, but I will have it once it finishes running again.

For reference I have 64GB ram.

Unclear above is the fact that vbebdot and omegabnbdot are HUGE equations (spanning several pages).