Your function, f = @(t,y) ... has vector y with just three components, but you call it in your Euler routine with i = 4 ...
y needs to be a matrix of 3xn where n is the number of points at which you want to evauate the function (which itself had a syntax error).
If you replace all the lines in your code, from that function definition, with the following code, your program will work. (You will, of course, have to replace the arbitrary constants/times/etc that I've used with your actual data).
f=@(t,y)[(F1*(Tc-y(1))+Fr*(T2_a-y(1))+(Q1/(rho_w*Cp)))/(V1);
(F1*(T1_a-y(2))+F2*(Tc-y(2))-Fr*(y(2)-T1_a)+(Q2-(2*pi*h2*R2*U*(y(2)-Ta))/rho_w*Cp))/(A2*h2);
(F1+F2-k*((Fd_A*((y(3)^3))+Fd_B*((y(3)^2))-Fd_C*(y(3))+0.10541)^0.5)*1e-03)/A2];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following data is arbitrary and needs to be replaced by true values
t0 = 0; tf = 1; % Initial and final times/positions respectively
n = 100; % Number of times/segments
h = (tf-t0)/n; % step size
t = t0:h:tf; % times
y = zeros(3,n+1); % Storage space
y(1,1) = Tc; y(2,1) = T1_a; y(3,1) = 1; % Initial values
% Simple Euler integration
for i = 1:n
y(:,i+1)=y(:,i)+h*f(t(i),y(:,i));
end
% Plot results
subplot(3,1,1)
plot(t,y(1,:)),grid
xlabel('t'),ylabel('y1')
subplot(3,1,2)
plot(t,y(2,:)),grid
xlabel('t'),ylabel('y2')
subplot(3,1,3)
plot(t,y(3,:)),grid
xlabel('t'),ylabel('y3')
