Speeding up automatic symbolic jacobian creation in Matlab

Question

I have a code that creates a set of equations based on nodal coordinates in cartesian space. I then create a jacobian based on these using automatic sym generation and a series of nested anonymous functions. I understand that this is not the most effective way to perform the operation, but my only other idea is to have .m files written and implemented at runtime. Here's a sample of the equations and methodology I use.

% c is created from a graph, here is the associated c for the example
c=[[1,-1,0,0,0,0];[0,1,-1,0,0,0];[0,0,1,-1,0,0];[0,0,0,1,-1,0];[0,0,0,0,1,-1];[1,0,0,0,0,-1];[0,1,0,0,-1,0];[1,0,0,0,-1,0];[0,1,0,0,0,-1];[0,1,0,-1,0,0];[0,0,1,0,0,-1]];
% base system, six nodes in 2d
xsys = [[0 1 2 2 1 0]' [1 1 1 0 0 0]'];
dexample = rand(6,2)*0.1;
% xsys is in this case an n by d matrix of coordinates (nodes by dimensions) 
%  and d is the displacement of xsys
u = @(d) c *(xsys(:,1)+d(:,1));
v = @(d) c *(xsys(:,2)+d(:,2));
U = @(d) diag(u(d));
V = @(d) diag(v(d));
Linv = @(d) diag(1./sqrt(U(d)*u(d)+V(d)*v(d)));
L0inv = Linv(xsys);
Q = @(d) EA*(L0inv - Linv(d))

% created for the sake of differentiation
dd = sym('d%d%d',size(xsys),'real');

% A and B here are sparse 0/1/-1 matrices. I am applying a load matrix
%  (same size as xsys) vector fe
eq1 = @(d) A' * B' * Q(d) * c * (xsys+d) + A * fe;

% eq1, in this particular case, creates a matrix of 8 equations using a
%  total of 12 variables contained in d
j1 = @(d) [reshape(diff(eq1(d),dd(1)),[8 1])';
           reshape(diff(eq1(d),dd(2)),[8 1])';
           reshape(diff(eq1(d),dd(3)),[8 1])';
           reshape(diff(eq1(d),dd(4)),[8 1])';
           reshape(diff(eq1(d),dd(5)),[8 1])';
           reshape(diff(eq1(d),dd(6)),[8 1])';
           reshape(diff(eq1(d),dd(7)),[8 1])';
           reshape(diff(eq1(d),dd(8)),[8 1])';
           reshape(diff(eq1(d),dd(9)),[8 1])';
           reshape(diff(eq1(d),dd(10)),[8 1])';
           reshape(diff(eq1(d),dd(11)),[8 1])';
           reshape(diff(eq1(d),dd(12)),[8 1])'];

% varargs2fun outputs argin as argout
objective = @(d) varargs2fun( eq1(d), j1(d) );

I then perform an optimization to set all values output from @eq1 to 0, based on the load application fe. This procedure takes 200 seconds, the majority of which is spent in mupadmex. Help? I need to drop the order of time by 5 or 6...

Again, my goal is for the jacobian to be automatically created, because the variables a b and c will change, and affect the number of equations in @eq1