The expression
[-ka3*A\[0; 1 ;0] zeros(3,1)]
must be 6x6 instead of 3x2 in order to work for your system.
And the numerical integration method is "Runge-Kutta", not "Range-Kutta". Two people, one with name "Runge", the other with name "Kutta".
clear all, close all, clc
s = 7.5; % semi span
c = 2; % chord
hline = 80; % hinge line as percentage of chord
xh = hline / 100 * c; % distance of hinge line aft of leading edge
fa = 40; % elastic axis as percentage of chord
xf = fa / 100 * c; % position of elastic axis aft of le
xac = 0.25 * c; % distance of ac aft of le
xe = xf - xac; % distance of elastic axis aft of ac
e = xe / c; % eccentricity of fa aft of ac per chord
heave = 3; % heave angle
alpha = 3; % pitch angle
beta = 3; % control surface angle
chat = 2 * (xh / c) - 1;
T4 = - acos(chat) + chat * sqrt(1 - chat ^ 2);
T10 = sqrt(1 - chat ^ 2) + acos(chat);
T11 = acos(chat) * (1 - 2 * chat) + sqrt(1 - chat ^ 2) * (2 - chat);
T12 = sqrt(1 - chat ^ 2) * (2 + chat) - acos(chat) * (2 * chat + 1);
% Aerodynamic coefficients
aw = 2 * pi; % lift per incidence
ac = aw * T10 / pi; % lift per control rotation
bw = e * aw; % pitching moment per control rotation
bc = e * aw * T10 / pi; % pitching moment per control rotation
cw = - T12 / 2; % hinge moment per incidence
cc = - T12 * T10 / 2 / pi; % hinge moment per control rotation
Mtdot = -1.2; % unsteady torsional damping term
Mbdot = -0.1; % unsteady control rotation damping term
rho = 1.225; % air density
% Mass data
m = 400; % mass per unit area
mmain = m; % mass per unit area of main surface
mcont = m; % mass per unit area of control surface
% Stiffness data
EI = 4e+7; % flexural rigidity
GJ = 8e+6; % torsional rigidity
kbeta = 1e+4; % control stiffness per unit span
kh = 4 * EI / s^3; % bending spring stiffness
ka = GJ / s; % torsion spring stiffness
ka3 = 100*ka*alpha^3; %non linear stiffness
GJstr = sprintf('%0.5g',GJ); EIstr = sprintf('%0.5g',EI);
astring = ['xf/c = ', num2str(xf/c),' xh/c = ',num2str(xh/c),' EI = ',EIstr,...
' GJ = ',GJstr,' kbeta = ', num2str(kbeta)];
disp(astring)
modes = 3;
% Inertia matrix
A = zeros(modes, modes);
% Main surface
Amain = zeros(modes, modes);
Amain(1,1) = xh * s / 5;
Amain(1,2) = (xh^2 / 2 - xh * xf) * s / 4;
Amain(2,1) = Amain(1,2);
Amain(2,2) = (xh^3 / 3 - xh^2 * xf + xf^2 * xh) * s / 3;
% Control surface
Acont = zeros(modes, modes);
Acont(1,1) = (c - xh) * s / 5;
Acont(1,2) = ((c^2 - xh^2) / 2 - xf * (c - xh)) * s / 4;
Acont(1,3) = ((c^2 - xh^2) / 2 - xh * (c - xh)) * s / 3;
Acont(2,2) = ((c^3 - xh^3) / 3 - xf * (c^2 - xh^2) + xf^2 * (c - xh)) * s / 3;
Acont(2,3) = ((c^3 - xh^3) / 3 - xf * (c^2 - xh^2) / 2 - xh * (c^2 - xh^2) / 2 + xf * xh * (c - xh)) * s / 2;
Acont(3,3) = ((c^3 - xh^3) / 3 - xh * (c^2 - xh^2) + xh^2 * (c - xh)) * s;
Acont(2,1) = Acont(1,2);
Acont(3,1) = Acont(1,3);
Acont(3,2) = Acont(2,3);
A = mmain * Amain + mcont * Acont;
% Structural damping matrix is zero
D = zeros(modes, modes);
B = zeros(modes, modes); % aero damping - based on rho*V*B
C = zeros(modes, modes); % aero stiffness - based on rho*V^2*C
B(1,1) = aw * c * s / 10;
B(1,2) = 0;
B(1,3) = 0;
B(2,1) = - bw * c^2 * s / 8;
B(2,2) = - Mtdot * c^3 * s / 24;
B(2,3) = 0;
B(3,1) = - cw * c^2 * s / 6;
B(3,2) = 0;
B(3,3) = - Mbdot * c^3 * s / 8;
C(1,1) = 0;
C(1,2) = aw * c * s / 8;
C(1,3) = ac * c * s / 6;
C(2,1) = 0;
C(2,2) = - bw * c^2 * s / 6;
C(2,3) = - bc * c^2 * s / 4;
C(3,1) = 0;
C(3,2) = - cw * c^2 * s / 4;
C(3,3) = - cc * c^2 * s / 2;
E = zeros(modes, modes);
E(1,1) = kh;
E(2,2) = ka;
E(3,3) = kbeta * s;
n = 3;
V = 117;
t = linspace(0,10,500);
Mat=[zeros(n,n),eye(n,n);
-A\(rho*V*V*C+E), -A\(rho*V*B+D)];
%X1 = [x(1)
% x(2)
% x(3)
% x(4)
% x(5)
% x(6)];
%x(1) = x1;
%x(2) = x2;
%x(3) = x3;
%x(4) = x4;
%x(5) = x5;
%x(6) = x6;
%não linearidade em alpha
%X2 = [0
% x(2)^3
% 0
% 0
% 0
% 0];
size(Mat)
size([-ka3*A\[0; 1 ;0] zeros(3,1)])
Mat3 = @(t,x) Mat*x + [-ka3*A\[0; 1 ;0] zeros(3,1)]*[0; x(2)^3;0;0;0;0];
x0 = zeros(6,1);
tspan = [0 1];
[T,X] = ode45(Mat3,tspan,x0)
plot(T,X)