Problem with Fplot : Array Input ?

Hello,

I have a problem with my Matlab Code, indeed fplot doesn't working and I have the following warning message :

 "Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output
with the same size and shape as the input arguments." 

 "The following error was reported evaluating the function in FunctionLine update: Matrix dimensions must agree."

I don't understand where I made a mistake... Should I use "ArrayFun" but I don't really know how does it works ? If someone can explain to me, it will be really nice !

Here is my code :

%Data
L_3D=4.19; %Ship length (m
Aw = 1; % wave amplitude (m)
T = 2*pi/7; % wave period(s)
%w = 2*pi/T ;% pulsation (rad/s)
L = 50; % Wavelength (m)
v= L/T; %Wave velocity
k = 2*pi/L; % L = 2*pi/k  L wavelength
rho = 1025 ; % density of seawater (kg/m^3)
%CM = 2 ; % inertia coefficient
%CD = 1 ; % drag coefficient
D1 =0.5; % cylinder diameter(m)
m = 2 ; % cylinder mass (kg)
S = pi*(D1/2)^2 ; %cylinder surface (m²)
g = 9.81 ; %gravity acceleration (m/s²)
y_a= 0.025;%Aw/2; Amplitude of heave (m)
x_a=L/2; %Amplitude of sway (m)
beta_a=1.5*pi/180; %Amplitude of roll(rad)
delta = A*L_3D*Sigma*0.68; %Displacement volume (m^3)
BG = 0.1875; %Vertical distance between centre of gravity and centre of buoyancy (m)*
GM = 0.17; %Metacentric height (m)
Fn = v/sqrt(g*L_3D); %Froude number

%Radiation potential

%For heave
tic
f2_h =@(w,x) -( sin(2*x) - ((w*w*M)/g)*( - sin(x) - (1/3)*a1*sin(3*x) + (3/5)*a3*sin(5*x))) + ( -(M/B)*((sin(-x) -a1*sin(x) + a3*sin(3*x)))*(-(w*w*M)/g)*(-1+(1/3)*a1+(3/5)*a3));
f4_h =@(w,x) -( sin(4*x) - ((w*w*M)/g)*( - (1/3)*sin(3*x) - (1/5)*a1*sin(5*x) + (3/7)*a3*sin(7*x))) + ( -(M/B)*((sin(-x) -a1*sin(x) + a3*sin(3*x)))*(-(w*w*M)/g)*((1/3)-(1/5)*a1-(3/7)*a3));

f22_h=@(w,x)(f2_h(w,x).*f2_h(w,x));
f24_h=@(w,x)(f2_h(w,x).*f4_h(w,x));
f44_h=@(w,x)(f4_h(w,x).*f4_h(w,x));

H=@(w,x) -(-(M/B)*((sin(-x) -a1*sin(x) + a3*sin(3*x))));

Phi_A02_h =@(w,x) ( cos(2*x) - ((w*w*M)/g)*( - cos(x) - (1/3)*a1*cos(3*x) + (3/5)*a3*cos(5*x)));
Phi_A04_h =@(w,x) ( cos(4*x) - ((w*w*M)/g)*( - (1/3)*cos(3*x) - (1/5)*a1*cos(5*x) + (3/7)*a3*cos(7*x)));
Ksi_A02_h =@(w,x) ( sin(2*x) - ((w*w*M)/g)*( - sin(x) - (1/3)*a1*sin(3*x) + (3/5)*a3*sin(5*x)));
Ksi_A04_h =@(w,x) ( sin(4*x) - ((w*w*M)/g)*( - (1/3)*sin(3*x) - (1/5)*a1*sin(5*x) + (3/7)*a3*sin(7*x)));
Phi_B0c_h= @(w,x) pi.*exp(-((w*w*D)/(2*g))*cos(x)).*cos(((w*w*D)/(2*g))*sin(x));
Phi_B0s1_h =@(w,x) pi.*exp(-((w*w*D)/(2*g))*cos(x)).*sin(((w*w*D)/(2*g))*sin(x));
Ksi_B0c_h= @(w,x) pi.*exp(-((w*w*D)/(2*g))*cos(x)).*sin(((w*w*D)/(2*g))*sin(x));
Ksi_B0s1_h= @(w,x) -pi.*exp(-((w*w*D)/(2*g))*cos(x)).*cos(((w*w*D)/(2*g))*sin(x));
Int_Phi_1_h= @(w,x) (integral(@(k) (((w*w)/g)*sin(k*(D/2).*cos(x)).*(exp(-k*(D/2).*sin(x)))-(k.*cos(k*(D/2).*cos(x))).*(exp(-k*(D/2).*sin(x))))./((k.*k)+((w*w)/g)*((w*w)/g)),0,Inf));
Phi_B0s_h=@(w,x) Phi_B0s1_h(w,x) + Int_Phi_1_h(w,x);
Int_Ksi_1_h= @(w,x) (integral(@(k) (((w*w)/g)*cos(k*(D/2).*cos(x)).*(exp(-k*(D/2).*sin(x)))+(k.*sin(k*(D/2).*cos(x))).*(exp(-k*(D/2).*sin(x))))./((k.*k)+((w*w)/g)*((w*w)/g)),0,Inf));
Ksi_B0s_h=@(w,x) Ksi_B0s1_h(w,x) + Int_Ksi_1_h(w,x);

I22_h =@(w) integral(@(x) f22_h(w,x),0,(pi/2));
I24_h =@(w) integral(@(x) f24_h(w,x),0,(pi/2));
I44_h =@(w) integral(@(x) f44_h(w,x),0,(pi/2));

g2_h=@(w,x) (Ksi_B0c_h(w,x)- H(w,x).*Ksi_B0c_h(w,pi/2));
g3_h=@(w,x) (g2_h(w,x).*f2_h(w,x));
g4_h=@(w,x) (g2_h(w,x).*f4_h(w,x));
g5_h=@(w,x) (Ksi_B0s_h(w,x)- H(w,x).*Ksi_B0s_h(w,pi/2));
g6_h=@(w,x) (g5_h(w,x).*f2_h(w,x));
g7_h=@(w,x) (g5_h(w,x).*f4_h(w,x));

%Calculation of the coefficients Pm and Qm

J1_h =@(w) integral(@(x) g3_h(w,x),0,(pi/2));
J2_h =@(w) integral(@(x) g4_h(w,x),0,(pi/2));
J3_h =@(w) integral(@(x) g6_h(w,x),0,(pi/2));
J4_h =@(w) integral(@(x) g7_h(w,x),0,(pi/2));

A1_h =@(w) [I22_h(w) I24_h(w); I24_h(w) I44_h(w)];
B1_h =@(w) [J1_h(w) ; J2_h(w)];
B2_h =@(w) [J3_h(w) ; J4_h(w)];

Pm_h=@(w)inv(A1_h(w))*B1_h(w);
Qm_h=@(w) inv(A1_h(w))*B2_h(w);

A10=[1,0];
A01=[0,1];

A0_h =(@(w) Ksi_B0c_h(w,pi/2) + (A10*Pm_h(w)).*Ksi_A02_h(w,pi/2) + (A01*Pm_h(w)).*Ksi_A04_h(w,pi/2));
A0_h(1)
B0_h =(@(w) Ksi_B0s_h(w,pi/2) + (A10*Qm_h(w)).*Ksi_A02_h(w,pi/2) + (A01*Qm_h(w)).*Ksi_A04_h(w,pi/2));

RAO=@(w) (pi*B.*w.*w)./(g.*sqrt((A0_h(w).*A0_h(w))+(B0_h(w).*B0_h(w))))
RAO(1)
fplot(RAO,[0,2*pi])
toc

Thanks,