How do I return a matrix in a function? I keep getting an error that says "Assignment has more non-singleton rhs dimensions than non-singleton subscripts"

Here is my code I wrote,

function [out] = two_link_bot(theta1, theta1_dot, theta1_ddot, beta, beta_dot, beta_ddot)
m1 = 2;
m2 = 1;
L1 = 2;
L2 = 1;
g = 9.81;
I1 = 0.667;
I2 = 0.083;
c1 = 0.005;
c2 = 0.005;
mee = 0.1;
mm2 = 0.25;
% Tf1 = 1.6;
% Tf2 = 0;
% x
% C1 = (m1*L1^2/4) + I1 + (m2*L1^2) + (mm2*L1^2) + (mee*L1^2);
% C2 = ((m2/2)+mee)*L1*L2;
% C3 = (m1*g*L1/2) + (m2*g*L1) + (mm2*g*L1) + (mee*g*L1);
% C4 = I2 + ((m2*L2^2)/4) + (mee*L2^2);
% C5 = (m2/2 + mee)*g*L2;
% C6 = (m1/2) + mm2;

% Will create a matrix A multiplied by matrix x that equals b, 
% X will contain our unkowns, so we can solve for x in the following order
% O1x O1y O2x O2y tau1 tau2

A = zeros(6);
B = zeros(6,1);

A(1,1) = 1;
A(1,3)= 1;
A(2,2) = 1;
A(2,4) = 1;
A(3,1) = (L1/2).*sin(theta1);
A(3,2) = -(L1/2).*cos(theta1);
A(3,3) = -(L1/2).*sin(theta1);
A(3,4) =(L1/2).*cos(theta1);
A(3,5) = 1;
A(3,6) = 1;
A(4,3) = -1;
A(5,4) = -1;
A(6,3) = -(L2/2)*sin(beta);
A(6,4) = (L2/2)*cos(beta);
A(6,6) = 1;
B(1) = -((m1/2) + mm2)*(L1*sin(theta1).*theta1_ddot) - ((m1/2) + mm2)*L1*cos(theta1).*theta1_dot.^2;
B(2) = C6(L1)*cos(theta1).*theta1_ddot + m1*g+mm2*g-C6*L1*sin(theta1).*theta1_dot.^2;
B(3) = (I1+mm2*L1^2/2)*theta1_ddot + c1*theta1_dot - c2*theta1_dot+mm2*g*L1/2*cos(theta1);
B(4) = -(m2+mee)*L1*sin(theta1).*theta1_ddot - (m2/2 + mee)*L2*sin(beta)-(m2+mee)*L1*cos(theta1)*theta1_dot.^2 - (m2/2+mee)*L2*cos(beta)*beta_dot.^2;
B(5) = (m2+mee)*L1*cos(theta1).*theta1_ddot + (m2/2+mee)*L2*cos(beta).*beta_ddot - (m2+mee)*L1*sin(theta1).*theta1_dot.^2 - (m2/2 +mee)*L2*sin(beta)*beta_dot.^2+m2*g+mee*g;
B(6) = mee*(L2/2)*L1*cos(beta-theta1).*theta1_ddot + (I2+mee*(L2^2)/2)*beta_dot + mee*g*(L2/2)*cos(beta) + c2*theta1_dot+mee*(L2/2)*L1*sin(beta-theta1).*theta1_dot.^2;

out = A\B;
end

There is another file that is a .mat file that has the values for the thetas and betas.

I've attached it.

I want to get a graph of tau1 and tau2 profiles as the angles are input.