Help solving system of triq equations

I am using matlab to create a script to solve for the angles theta 3 and theta 4 with an input of theta 1 and theta 2 of a fourbar

theta 1 is set to 0 for this data set however next time i call this function it will not be zero and will refrence the output from the first iteration

%Input Linkage lengths
d=11.875;
a=2.5;
b=11.875;
c=3.75;
r5=5;
r6=6;
r7=7;
r8=8;

L1Theta = zeros(30,4);
L1Theta(:,2) = [linspace(0,2*pi,length(L1Theta))];


positionsolver(a,b,c,d,L1Theta(:,1),L1Theta(:,2),length(L1Theta))


    function [theta3,theta4] = positionsolver(A,B,C,D,theta1,theta2,L)
        for n = 1:L
            Theta1 = theta1(n)
            Theta2 = theta2(n)
            syms Theta1 Theta2 A B C D Theta3 Theta4

            eq1 = (A*cos(Theta2(n)) + B*cos(Theta3) - C*cos(Theta4) - D*cos(Theta1(n))) == 0
            eq2 = (A*sin(Theta2(n)) + B*sin(Theta3) - C*sin(Theta4) - D*sin(Theta1(n))) == 0
            vars = [Theta3,Theta4]
            sol = solve([eq1,eq2],vars)
            theta3(n) = simplify(sol.Theta3)
            theta4(n) = simplify(sol.Theta4)
        end

    end


end

Answer 1

[Torsten]

%Input Linkage lengths
d=11.875;
a=2.5;
b=11.875;
c=3.75;
r5=5;
r6=6;
r7=7;
r8=8;
L1Theta = zeros(30,4);
L1Theta(:,2) = [linspace(0,2*pi,length(L1Theta))];

theta1 = L1Theta(:,1);
theta2 = L1Theta(:,2);
L = length(L1Theta);

syms A1 A2 B C real
syms Theta3 Theta4
eq1 = B*cos(Theta3) - C*cos(Theta4) == A1;
eq2 = B*sin(Theta3) - C*sin(Theta4) == A2;
[Theta3,Theta4] = solve([eq1,eq2],[Theta3 Theta4])
Theta3 = simplify(Theta3)
Theta4 = simplify(Theta4)

% Compute results
for n = 1:L
   a1 = d*cos(theta1(n))-a*cos(theta2(n));
   a2 = d*sin(theta1(n))-a*sin(theta2(n));
   theta3(n,:) = double(subs(Theta3,[B,C,A1,A2], [b,c,a1,a2]));
   theta4(n,:) = double(subs(Theta4,[B,C,A1,A2], [b,c,a1,a2]));
end
theta3
theta4

% Check results
for n = 1:L
  a1 = d*cos(theta1(n))-a*cos(theta2(n));
  a2 = d*sin(theta1(n))-a*sin(theta2(n));
  res1(n,:) = b*cos(theta3(n,:))-c*cos(theta4(n,:))-a1;
  res2(n,:) = b*sin(theta3(n,:))-c*sin(theta4(n,:))-a2;
end
res1
res2