I'm not sure how the equation is supposed to be defined but the main problem here is that you're multiplying arrays of different dimensions. If instead of generating 10 random numbers you generate 51 (and make a little adjustment in t_tank3) , the code will work. Maybe it worked before with only 2 dimensions because matlab would expand the array properly, while now you have conflicting dimensions.
L=400; %km
v_t = unifrnd(30,50,51,1); % 10 to 51 for matching dimensions
v_tr = unifrnd(90,110,51,1);
v_tt = unifrnd(70,90,51,1);
t_l= unifrnd(0,(1/6),51,1);
t_ul= unifrnd(0,(1/6),51,1);
%
% compute equations for multiple x,y values
[x,y]=meshgrid(30:50,100:150);
t_tank1=(t_l+x./v_tt+t_ul)+(L-x)./v_t;
t_tank2=(t_l+x./v_tt+t_ul)+x./v_tr +(t_l+y./v_tt+t_ul)+(L-y)./v_t;
%here L./v_tt, to make it element-wise
t_tank3=(t_l+x./v_tt+t_ul)+x./v_tr +(t_l+y./v_tt+t_ul)+y./v_tr+(t_l+L./v_tt+t_ul);
t_operationpart=max(t_tank1,t_tank2);
t_operation=max(t_operationpart,t_tank3);
% plot t_operation
surf(x,y,t_operation)
xlabel('x (km)')
ylabel('y (km)')
